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Abstract 


Pulsar timing arrays (PTAs) are galactic-scale gravitational wave (GW) detectors. Each individual arm, composed 
of a millisecond pulsar, a radio telescope, and a kiloparsecs-long path, differs in its properties but, in aggregate, can 
be used to extract low-frequency GW signals. We present a noise and sensitivity analysis to accompany the 
NANOGrav 15 yr data release and associated papers, along with an in-depth introduction to PTA noise models. As 
a first step in our analysis, we characterize each individual pulsar data set with three types of white-noise 
parameters and two red-noise parameters. These parameters, along with the timing model and, particularly, a 
piecewise-constant model for the time-variable dispersion measure, determine the sensitivity curve over the low- 
frequency GW band we are searching. We tabulate information for all of the pulsars in this data release and present 
some representative sensitivity curves. We then combine the individual pulsar sensitivities using a signal-to-noise 
ratio statistic to calculate the global sensitivity of the PTA to a stochastic background of GWs, obtaining a 
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minimum noise characteristic strain of 7 x 107 


at 5 nHz. A power-law-integrated analysis shows rough 


agreement with the amplitudes recovered in NANOGrav's 15 yr GW background analysis. While our 
phenomenological noise model does not model all known physical effects explicitly, it provides an accurate 
characterization of the noise in the data while preserving sensitivity to multiple classes of GW signals. 


Unified Astronomy Thesaurus concepts: Millisecond pulsars (1062); Pulsar timing method (1305); Gravitational 
wave astronomy (675); Gravitational wave detectors (676); Radio astronomy (1338) 


1. Introduction 


Millisecond pulsars (MSPs) are extremely stable rotators, 
with long-term stability comparable to that of accurate atomic 
clocks (Matsakis et al. 1997; Hobbs et al. 2012, 20202). They 
are therefore uniquely well-suited astronomical objects to probe 
the low-frequency window of gravitational waves (GWs), from 
nanohertz to microhertz (Sazhin 1978; Detweiler 1979; Hel- 
lings & Downs 1983; Foster & Backer 1990; Perera et al. 
2018). 

The most promising sources of low-frequency GWs are 
inspiraling supermassive black hole binaries (SMBHBs). The 
superposition of the GWs produced by the SMBHB population 
forms the stochastic GW background (GWB), and the most 
massive, closest binaries could be resolved as continuous GWs 
produced by individual sources (see Taylor 2021, for review). 
The North American Nanohertz Observatory for Gravitational 
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Waves (NANOGrav; Alam et al. 2021a) Collaboration 
regularly monitors an array of highly stable MSPs spread 
across the sky to achieve the submicrosecond timing precision 
required for GW detection. In a GWB search, a common long- 
time-correlated (with more power at low frequency, or “red’’) 
signal should be detected among the pulsars in the array, 
combined with a quasi-quadrupolar signature in the angular 
correlations between pulsar pairs predicted by general 
relativity, the so-called Hellings—Downs correlations (Hellings 
& Downs 1983). NANOGrav reported the detection of a 
common red signal in our 12.5 yr pulsar timing data set; 
however, the observations were insufficient to exhibit the 
Hellings—Downs correlations required to definitively associate 
this signature with GWs (Arzoumanian et al. 2020). Similar 
GW-search results have subsequently been reported by other 
pulsar timing array (PTA) collaborations, including the 
European PTA (EPTA; Chen et al. 2021), the Parkes PTA 
(PPTA; Goncharov et al. 2021b), the International PTA (IPTA; 
Antoniadis et al. 2022; a collaboration between NANOGrav, 
the EPTA, the PPTA), and the recently formed Indian PTA 
(InPTA; Tarafdar et al. 2022). 

In this series of papers presenting the NANOGrav 15 yr data 
release, we report statistically significant evidence for Hellings 
—Downs spatial correlations in the timing residuals for an 
ensemble of 67 pulsars (see Agazie et al. 2023b, 
hereafter NGl15gwb). Because of the stochastic nature of the 
signal, a careful noise analysis of our experiment is imperative. 
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As we have done in our previous papers (Arzoumanian et al. 
2016, 2018, 2020), we explicitly model both white (time- 
uncorrelated) and nonwhite (time-correlated) noise components 
present in the detector and incorporate those results in 
additional analysis steps. Unlike in our previous data releases, 
in which this topic has been discussed in either the data 
presentation or GWB search papers, we devote a separate paper 
to the noise analysis of our 15 yr data release (see Agazie et al. 
2023a, hereafter NG15). This paper describes and employs a 
noise analysis technique specifically designed for the data in 
this release. Our approach has been developed over our last 
four data releases and incorporates empirical expressions of 
noise sources discussed in Cordes & Shannon (2010), Cordes 
(2013), Stinebring (2013), Demorest et al. (2013), Arzouma- 
nian et al. (2015), Arzoumanian et al. (2016), Lam et al. (2016), 
and Lam et al. (2017). Furthermore, we explicitly outline the 
structure of the covariance matrix we use since it incorporates 
our knowledge of noise source variances and their relation, if 
any, to each other. 

The structure of the paper is as follows. In Section 2 we 
present the known astrophysical sources of noise in our 
experiment. We present in Section 3 a phenomenological noise 
model that characterizes the white noise (WN) and its 
covariance, as well as developing a simple spectral model for 
red noise (RN) present in the data, which includes possible 
contributions from the GWB and other sources. The covariance 
matrices that take this knowledge of the noise and propagate it 
into later analyses are described in Section 4. Noise 
characterization and its effect on our sensitivity to a GWB 
are summarized through transmission functions and sensitivity 
curves in Section 5. Finally, we consider possible future 
improvements to noise characterization and reduction in 
Section 6. Notation used throughout the paper is listed in 
AppendixB. 


2. Astrophysical Impacts on Pulse Arrival Times 


Pulsar timing at the submicrosecond level relies on using 
MSPs with extremely stable rotation rates as high-precision 
cosmic clocks. We measure times of arrival (TOAs) of pulses 
from a pulsar by producing a summed pulse profile at every 
observation epoch and convolving with a pulse profile template 
using a matched-filter approach where the observed profile is a 
scaled and shifted version of a perfect template shape with 
additive noise (e.g., Taylor 1992). In this way, we are able to 
track a reference longitude on the pulsar, an arbitrarily chosen 
reference meridian. Measured TOAs are then subtracted from 
those predicted by a timing model to calculate pulsar timing 
residuals. Many known sources of noise both intrinsic and 
extrinsic to the pulsar impact the measured TOAs and must be 
carefully modeled in order to detect the timing signatures 
expected as a result of GWs in the residuals. 

In order to keep this paper self-contained, we briefly discuss 
the numerous astrophysical perturbations to the TOAs, the 
fundamental measured quantities in our experiment. This is a 
deep, rich, and complicated subject. See Cordes & Shannon 
(2010), Cordes (2013), Stinebring (2013), and Verbiest & 
Shaifullah (2018) for more detail on all of these effects. Each 
subsection is followed by one or a few labels, such as 
“chromatic,” “achromatic,” etc., describing the noise. Chro- 
matic (or chromaticity) refers to a dependence on the radio 
frequency of the observations. Achromatic features lack this 
dependence. Spatially correlated noise processes show 
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statistical cross-correlations that may depend on the sky 
position of the pulsars or possess the same correlation across 
all pulsars. These noise processes are particularly important to 
understand since PTAs search for the stochastic GWB, which is 
expected to be spatially correlated according to the Hellings 
—Downs relationship. 

In the following subsections noise sources are discussed 
roughly in order from emission mechanisms at the pulsar, 
following the propagation order, through to observation of 
radio pulses at the telescopes. See Figure 1 for a visual 
summary of these effects and their physical location relative to 
the pulsar— Earth lines of sight. 


2.1. TOA Errors from the Pulsar Emission Process, Including 
Jitter—Varying Chromaticity 


Phenomenologically, pulse jitter comprises the stochastic 
phase variation of single pulses with respect to a mean phase 
that is locked to the rotation of a neutron star. It appears to be a 
ubiquitous feature of radio pulsars and magnetars but is a much 
smaller (f not absent) effect in high-energy emission of 
pulsars. Pulse jitter limits the fundamental precision of our 
TOAs at the emission region itself. While the rotation of the 
neutron star represents the clock mechanism for our experi- 
ment, the TOAs we measure correspond to radio emission 
generated in the pulsar's magnetosphere. Pulse shapes averaged 
over many rotations are generally very stable, allowing for 
TOA measurements to precisions O(10-?) the duration of pulse 
widths. However, the shapes of individual pulses are not stable. 
This long-studied phenomenon is known as pulse jitter, in 
which single-pulse components can vary in phase and 
amplitude, with phase variations of order the pulse width 
(Cordes & Downs 1985; Lam et al. 2016; Parthasarathy et al. 
2021). The net effect is that any finite sample of pulses will 
slightly differ in shape compared to another sample, creating a 
source of uncertainty in the measured phase of the spot of 
longitude that should be uncorrelated between observing 
epochs and scale with the number of observed pulses NW, as 
xl / VN . Jitter is strongly correlated across radio frequencies, 
though since pulse profiles vary in shape across frequency, the 
combined effect on timing due to jitter can have some 
frequency dependence, which has been measured for a large 
fraction of our MSPs (Lam et al. 2019; Hebel 2022). However, 
jitter is entirely uncorrelated from pulse to pulse, i.e., it acts as 
WN in the residuals (Cordes & Downs 1985). 

Giant radio pulses (GRPs) are individual pulses that have 
been observed to occur not periodically and that show 
amplitudes many times larger than average and amplitude 
distributions that follow a power law. They are another 
phenomenon that arises in a small subset of known pulsars. 
The MSP B1937+21 was the second emitter of GRPs detected 
after the Crab pulsar (Staelin & Reifenstein 1968; Wolszczan 
et al. 1984). GRPs occur in specific ranges of pulse phase and 
are thus locked to the spin phase of the neutron star, but their 
large amplitudes occur stochastically without any distinct 
temporal pattern. In the case of B19374-21, its GRPs have 
pulse widths ranging from microseconds to nanoseconds, with 
intensities that follow a power-law distribution and show 
brightness temperatures ranging up to 10? K for a GRP of 
15ns pulse width (Sallmen 1995; Cognard et al. 1996; 
Kinkhabwala & Thorsett 2000; Soglasnov et al. 2004). GRPs 
from B1937--21 have been observed to have little effect on the 
timing precision of this pulsar. Extensive studies performed by 
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Figure 1. An illustration of the signal path from the pulsar to the data products, highlighting a few relevant sources of noise. The pulsar emission itself is subject to 
jitter (and, in rare cases, glitches), and the source model includes fits for spin and spin-down, astrometric parameters, and possible binary orbit models. Propagation 
through the ISM imposes pulse dispersion and scattering, which are both frequency and time dependent. At the receiver, along with thermal noise, there are effects 
related to clock corrections and solar system barycenter parameters, among others, that need to be modeled in order to produce a data set that includes pulse TOAs and 


time-dependent DM estimates. 


McKee et al. (2019) show that the GRPs from B1937+21 
provide a timing solution negligibly different from the general 
timing solution of the pulsar. Exhaustive studies of lower-level 
GRP emission have not been performed for the remaining 
NANOGrav pulsars, though examination of several individual 
sources is underway. A frequently discussed pulsar character- 
istic for the occurrence of GRPs is a high magnetic field at the 
light cylinder (Cognard et al. 1996). Currently, it is not clear 
whether this characteristic is connected to the physics of GRPs. 
Further investigations need to be carried out before corresp- 
onding interpretations can be made. Another pulsar that also 
displays a high magnetic field at its light cylinder and was 
reported to be the third GRP emitter is PSR B1821—24A 
(Romani & Johnston 2001). Albeit not observed by NANO- 
Grav, it is observed by other PTAs. No studies of the effect of 
its GRPs on its timing precision have been published yet. 
Although caused by propagation effects in circumpulsar 
environments rather than the emission process itself, highly 
magnified pulses were discovered in the MSP B19574-20 
(Main et al. 2018) and have been observed in a variety of 
eclipsing MSPs (see Lin et al. 2023). Their effect, if any, on 
MSP timing precision is largely unexplored. 

Mode changing is a phenomenon in which a pulsar's average 
profile switches between two or more discrete modes over 
a variety of sometimes-periodic timescales, ranging from 


seconds to days. The most dramatic, and first identified, 

example of mode changing is the abrupt cessation of a pulsar's 

radio emission for several consecutive pulse periods, known as 

nulling and first identified by Backer (1970). Nulling has been 

observed in a large number of predominantly slow-period pulsars 

(Huguenin et al. 1970; Page 1973; Hesse & Wielebinski 1974). 

Nulling has not been observed in any MSP, but it is difficult 
to detect owing to the high sensitivity required to detect MSP 
single pulses. However, mode changing has been observed in 
three MSPs so far—PSR B19574-20 (Mahajan et al. 2018), 

PSR J0621--1002 (Wang et al. 2021), and PSR J1909—3744 
(Miles et al. 2022). Of the three, only PSR J1909—3744 is 
observed by NANOGrav. The mode changes detected from 
PSR J1909-3744 are characterized by a weak mode and a 
strong mode that are offset in phase and occur on a single- 
pulse timescale (Miles et al. 2022). Additionally, PSR J1713 
+0747 shows evidence of distinct modes of low-amplitude 
drifting subpulses (Liu et al. 2016). For both pulsars, these 
variations should contribute similarly to WN as pulse-to-pulse 
jitter. We have shown that noise on timescales of several 
hours was shown to integrate down as WN (Dolch et al. 2014; 
Shapiro-Albert et al. 2020). Mode changing would likely be 
detectable in a larger sample of the MSPs we observe with 
sufficient sensitivity. A detailed investigation of the impact of 
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the nulling/mode changing over the time spans we observe is 
ongoing. 


2.2. Rotational Irregularities (Spin Noise)—Achromatic 


Spin noise results from rotational irregularities in the pulsars. 
The interiors of MSPs are thought to consist of a differentially 
rotating superfluid core surrounded by an iron-rich crust 
(Langlois et al. 1998). Rotational irregularities are thought to 
arise from torques in the pulsar magnetosphere (Cheng 1987; 
Kramer et al. 2006; Cordes & Shannon 2008; Lyne et al. 2010; 
Gao et al. 2016) or coupling between the crust and superfluid 
core (Jones 1990; Melatos & Link 2014; van Eysden & 
Link 2018). The result is a stochastic variation in the pulse 
period detectable over timescales of many years. Fortunately, 
spin noise is positively correlated with spin frequency 


derivative, with an rms contribution of opin e over long 
timescales (Shannon & Cordes 2010; Lam et al. 2017; Lower 
et al. 2020). Since MSPs have very small frequency 
derivatives, they typically show negligible spin noise. One 
counterexample is the NANOGrav pulsar B1937--21 (Shannon 
& Cordes 2010). Beyond rotational irregularities from neutron- 
star-specific origins, another proposed cause of the observed 
timing variations of this pulsar is an asteroid belt of mass 
<0.05 Mẹ (Shannon et al. 2013). 

Glitches represent a special case of spin noise where the spin 
of a pulsar experiences a step change (see, e.g., Espinoza et al. 
2011). Although this usually results in a very small increase in 
spin frequency, examples of large increases (“slow glitches"; 
Shabanova 2005) have been reported, or even decreases (so- 
called “anti-glitches”) in X-ray pulsars and magnetars (Archi- 
bald et al. 2013; Ray et al. 2019). The glitch rate and typical 
fractional change of spin frequency show strong inverse 
correlations with pulsar age (Fuentes et al. 2017). They are 
therefore rare in MSPs, although two examples have been 
reported: first in the globular cluster pulsar B1821—24A 
(Cognard & Backer 2004), and later in the NANOGrav pulsar 
J0613—0200 (McKee et al. 2016). The fractional changes in 
spin frequency for these two glitches are only (8 + 1) x 10? 
and (2.521) x 10^, respectively, making these the two 
smallest glitches listed in the Jodrell Bank Glitch Catalogue®® 
(Espinoza et al. 2011). The low glitch amplitudes make glitches 
very hard to detect in MSPs, with the PSR J0613—0200 glitch 
occurring in early 1998 but not being reported until 2016. Note 
that this glitch occurred prior to the start of our data set. 
However, other small and unrecognized glitches from this or 
other MSPs could prove problematic for detection of some 
proposed GW targets of PTAs, as a GW memory burst 
(affecting only the pulsar term) is predicted to induce an 
apparent spin frequency step change in the timing residuals of a 
single pulsar that would be indistinguishable from the signature 
of a small glitch (Cordes & Jenet 2012). 


2.3. Profile Changes in Frequency and Time—Chromatic 


Pulse shapes intrinsically vary as a function of frequency 
and, for the most distant pulsars, as a result of multipath 
propagation through the interstellar medium (ISM), or scatter- 
ing; see Section 2.5. We account for this effect in two different 
ways. Our traditional narrowband timing method uses a single 
template to fit pulses in a number of discrete frequency 


68 https: / [www jb.man.ac.uk/pulsar/glitches.html 


Agazie et al. 


channels across a given frequency band and then corrects for 
constant time offsets between channels in the timing models 
with a fitted functional form of the offsets versus frequency 
(Agazie et al. 2023a). Our “wideband” method (Pennucci et al. 
2014; Pennucci 2019) uses a pulse portrait that contains 
information about the frequency-dependent pulse shape to fit 
for the TOA and the time-dependent dispersion measure (DM; 
the integrated column density of free electrons along the line of 
sight) delay simultaneously per epoch. 

The TOA creation process assumes that the pulse profile of 
the pulsar is constant and reproducible when integrating over a 
suitably long time (at least thousands of pulses). However, 
subtle pulse shape variations, likely attributable to either 
propagation effects or incorrect polarization calibration, have 
been observed for several NANOGrav pulsars over long 
timescales (Brook et al. 2018). 

In addition, in early 2021 (after the time span of our data 
set), the NANOGrav pulsar J1713+0747 was found to have 
experienced a drastic change in its pulse shape on a timescale 
of less than 1 day (Singha et al. 2021a; Lam 2021; Meyers & 
Chime/Pulsar Collaboration 2021; Xu et al. 2021), before 
relaxing back toward its original pulse shape over the course of 
several months (Jennings et al. 2022b). As the signature in the 
timing residuals (generated using templates derived from pre- 
event data) was similar to that expected from ISM propagation 
effects (Lin et al. 2021), this shape change was initially 
misinterpreted as being accompanied by a step change in DM. 
It was later demonstrated that, as the pulse shape evolution over 
frequency was different before and after the event (Singha et al. 
2021b), it is not trivial to measure a change in DM between 
epochs (Jennings et al. 20222) and that the observed chromatic 
behavior is entirely consistent with a frequency-dependent 
pulse shape change unaccompanied by changes in ISM 
propagation. This event has necessitated reexamination of the 
previous chromatic timing features in PSR J1713--0747 data 
noted by Lam (2021), as well as one in PSR J1643— 1224 
(Shannon et al. 2016); it is not yet clear how much of an impact 
this phenomenon has on our ability to model pulsar timing 
behavior. 


2.4. Orbital Irregularities—Achromatic 


The majority of pulsars observed by NANOGrav are in 
binary systems, typically with a white dwarf companion 
(Fonseca et al. 2016; Agazie et al. 2023a). A majority are in 
orbits well described by Keplerian parameters, as well as post- 
Keplerian parameters in many cases. However, several pulsars 
require different modeling of their binary systems, which we 
comment on here. Four (PSR J0023--0923, PSR J06364-5128, 
PSR J1705—1903, and PSR J22144-3000) are in low- 
eccentricity black widow systems where tidal and wind effects 
can cause measurable orbital variations. Three of these pulsars 
(excluding PSR J2214--3000) are modeled by higher-order 
orbital frequency derivatives. While PSR J1705—1903 is newly 
added to the NANOGrav PTA and several new higher-order 
orbital frequency derivatives have been measured in NG15 for 
the others, the orbital parameters for the other three measured 
in common with our 12.5 yr data set (Alam et al. 2021a) remain 
stable with additional data (Agazie et al. 2023a), suggesting 
that excess noise from parameter misestimation is low. Only 
PSR J1705—1903 shows significant RN, with a very shallow 
spectral index, suggesting that mismodeling does not signifi- 
cantly affect our GW analyses. Note that NANOGrav selects 
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for “well-behaved” systems from a timing perspective, and so 
other such black widow systems may indeed show significantly 
more noise owing to irregular angular momentum transfer, 
although since the orbital periods are all short compared to our 
GW signal (<Iday vs. years) only a small amount of 
sensitivity to GWs is lost overall (Bochenek et al. 2015). 

Beyond these systems, PSR J1024—0719 is a pulsar in a 
wide-binary system (Bassa et al. 2016; Kaplan et al. 2016) such 
that we cannot feasibly measure a complete orbit; instead, we 
model the orbital motion as a second derivative in the pulsar's 
spin frequency. In this case, the pulsar shows significant RN in 
individual noise and common noise analyses (see Table 2), and 
so unmodeled orbital variations may still be contributing to the 
noise in this pulsar. Lastly, there are two known pulsars in 
triple systems (Thorsett et al. 1999; Ransom et al. 2014), one of 
which (PSR J0337--1715) is currently being timed by 
NANOGrav owing to its high-precision TOAs but is not 
included in NG15. The timing of this pulsar requires higher- 
order effects from general relativity to be included in its timing 
model, and so the modeling procedure involves more costly 
likelihood evaluations. Nonetheless, even though it is a 
complicated system, only a small amount of excess noise has 
been measured beyond the template-fitting uncertainties 
(Archibald et al. 2018), with no RN and a simplified DM fit 
compared to our procedure (see next subsection), making triple 
systems potentially significant contributors to future PTA 
efforts. 


2.5. ISM Propagation Effects—Chromatic 


Dispersion is the dominant frequency-dependent propaga- 
tion effect in pulsar timing resulting from radio pulses traveling 
through the ionized ISM, interplanetary medium (IPM; or solar 
wind), and even Earth's ionosphere, with emission being 
temporally delayed as a function of frequency v by an amount 
proportional to DM/ 1°. The relative motions of Earth, the 
pulsar, and the ISM cause the sampled free electron content to 
vary, requiring us to estimate a DM for every observing epoch. 
These dispersive delays are also covariant with other propaga- 
tion effects and the frequency dependence of the pulse shape 
itself. While dispersive delays are accounted for in our timing 
models, other time-dependent propagation effects are not. 

Multipath propagation through the ISM causes extra time 
delays that may vary differently from the v ? dispersive delay 
along with distortions to the pulse shape that can cause 
misestimation of the TOAs in a variety of ways. The most 
visible effect is scattering, manifesting as a broadening of the 
pulse shape mathematically described as a convolution with a 
pulse broadening function (PBF). The PBF has an approxi- 
mately exponential tail, with the broadening scaling approxi- 
mately as » ^^, though the value of the index is highly 
dependent on the physics and geometry of the medium. In the 
low-scattering regime, the primary impact of scattering is to 
delay the pulse with the roughly same frequency scaling 
(Hemberger & Stinebring 2008). However, several of the 
NANOGrav pulsars show prominent scattering tails, especially 
those at high DM and/or those observed at our lowest radio 
frequencies. 

Another frequency-dependent effect results from the ray 
paths at different radio frequencies traversing different sets of 
electrons, resulting in slightly different DMs at different 
frequencies (Cordes et al. 2016). At sufficiently high levels 
of precision, this will require a  frequency-dependent 
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characterization of the DM at each epoch, currently outside 
the scope of our analysis. 

Refraction also results from multipath propagation and the 
deflection of the bulk set of rays through the ISM, causing the 
observed sky position of the pulsar to vary. This results in an 
additional geometric delay to the pulses proportional to v 4, 
with an additional delay caused by the variation in the angle of 
arrival proportional to v ? (Foster & Cordes 1990), the latter 
entirely covariant with the dispersive delay. 

Scintillation arises because pulsar images are extremely 
angularly compact, and the multiple ray paths interfere 
coherently with one another. The effect is most easily 
visualized observationally as heavy modulation of the dynamic 
spectrum, or the intensity of the pulsar as a function of 
frequency and time. Normally referred to as diffractive 
interstellar scintillation (DISS), the process at a particular 
epoch can be characterized by a coherence bandwidth, Avy, 
and a coherence timescale, At,, with both quantities inversely 
proportional to pulsar distance. DISS changes the PBF 
stochastically, resulting in a limit to the fundamental precision 
of our TOAs from propagation effects (vs. intrinsically at the 
neutron star for jitter), with larger uncertainties at lower 
frequencies (see, e.g., in Lam et al. 2018b). When the total 
observation time, Tobs, and total bandwidth, B, are very large 
with respect to At; and Avy, respectively, the effect of DISS on 
the timing is reduced. However, most NANOGrav MSPs are 
nearby, leading to coherence timescales of the same order or 
longer than our approximately 30-minute observation lengths 
and coherence bandwidths only factors of a few smaller than 
total observing bandwidths. This results in a “finite scintle” 
effect. This is independent of pulse signal-to-noise ratio (S/N) 
to first order (Cordes et al. 1990; Lam et al. 2016), making the 
rms timing noise due to scintillation covariant with jitter in that 
respect. However, its strong dependence on frequency allows it 
to be disentangled. In addition, measurements of the scintilla- 
tion parameters can be used as priors on the analysis, though 
these are not applied to our current analyses. 


2.6. Solar Wind Effects—Chromatic (Spatially Correlated) 


The solar wind is a stream of charged particles that escape 
from the solar corona owing to their high (~1 keV) energies 
(Marsch 2006). This leads to an ionized IPM with an electron 
number density that decays with increasing distance from the 
Sun. Throughout Earth’s orbit the line of sight to a given pulsar 
cuts through different parts of the IPM, resulting in an annual 
contribution to the pulsar's DM. The DM variations caused by 
the solar wind are subsumed in the generic DM variation 
model, DMX, used by NANOGrav; however, this excess DM 
can be modeled by assuming a spherically symmetric electron 
number density as a function of Sun— pulsar separation angle 0 
(Edwards et al. 2006), 


DM,,(0) = (4.85 x 1079 pc)no (1) 


sinQ 
where ng is the nominal electron number density at a distance 
of lau, typically assumed to be several particles per cm? in 
pulsar timing codes and measured to be in the same range (e.g., 
Madison et al. 2019; Hazboun et al. 2022). The maximum 
amplitude (i.e., at the smallest angular separation from the Sun) 
of this periodic contribution to the DM depends on the pulsar's 
angle from the ecliptic plane, with pulsars very close to the 
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ecliptic (latitude |8| 10?) showing very sharp peaks around 
the minimum solar separation (Jones et al. 2017; Donner et al. 
2020). The solar wind therefore represents both a chromatic 
and spatially correlated signal (Tiburzi et al. 2016) among 
pulsars in a timing array. 

However, the solar wind is neither spherically symmetric nor 
static over the course of the 11 yr solar cycle (e.g., Issautier 
et al. 2001), making the simple spherical model in Equation (1) 
inadequate to describe the solar wind contribution to DM for 
pulsars at our DM precisions (Madison et al. 2019; Tiburzi 
et al. 2019; Hazboun et al. 2022). The model also does not 
differentiate between the densities of the "fast" and "slow" 
solar wind (Tiburzi et al. 2021), leading to large annual 
changes in amplitude. As shown in Hazboun et al. (2022), our 
current practice of measuring and removing the DM delay on a 
per-observation basis allows us to account for the changing 
DMs along each line of sight except in extreme circumstances 
where the Sun— pulsar separation angle is extremely small and 
dual-receiver observations are not available. NANOGrav 
removes data for which the DM variation is likely to be too 
rapid to reliably represent it with DMX segments (see Section 
4.1 in NGI5 for more details), but it retains the raw data from 
observing pulsars near the cusp of closest line-of-sight 
approach to the Sun for solar wind studies. 


2.7. Solar System Ephemeris—Achromatic (Spatially 
Correlated) 


TOAs are measured at an observatory and then must be 
transformed to the quasi-inertial reference frame of the solar 
system barycenter (SSB; see, e.g., Lorimer & Kramer 2005 for 
full details). This is done through the use of a planetary 
ephemeris, notably the Development Ephemeris (DE; Stand- 
ish 1982) maintained by NASA’s Jet Propulsion Laboratory, 
most recently DE441 (Park et al. 2021), and the Observatoire 
de Paris—maintained INPOP ephemerides (Fienga et al. 2009), 
most recently INPOP19a (Fienga et al. 2019). The dominant 
source of uncertainty in these ephemerides arises from 
inaccurate measurements of planetary masses, particularly the 
outer planets, each of which contributes a sinusoid to the timing 
residuals with period equal to the planetary orbit. This 
signature will be spatially correlated among pulsars, roughly 
as a dipolar signal, in a PTA (Tiburzi et al. 2016). The orbital 
periods of the giant planets range from 11.9 to 164.8 yr, 
corresponding to frequencies of ~0.4—2.7 nHz, similar to the 
GWB frequency range that is probed by PTAs. The impact on 
PTA analyses is modeled by perturbing the orbital parameters 
in the ephemeris used to derive PTA noise limits and 
comparing with other models, as detailed in, e.g., Vallisneri 
et al. (2020) and Chen et al. (2021). The dependence of the GW 
Statistics in NANOGrav’s 11 yr data set on the version of the 
ephemeris used was the catalyst for the development of these 
methods (Arzoumanian et al. 2018). These models allow the 
analyses to be bridged between different versions of the solar 
system ephemeris to obtain equivalent results. Arzoumanian 
et al. (2020) and Vallisneri et al. (2020) showed a diminishing 
dependency on the choice of ephemeris for PTA GWB results, 
and all tests of different ephemeris versions on NGI5, 
including using a Bayesian ephemeris model (Vallisneri et al. 
2020), have shown insignificant differences in parameter 
recovery. In addition, PTA data sets have been used to identify 
non-GW common noise signals and have enabled limits to be 
placed on the masses of the outer planets, large asteroids, and 


Agazie et al. 


undiscovered planets in the outer solar system (Champion et al. 
2010; Caballero et al. 2018; Guo et al. 2019). 


2.8. Clock Errors—Achromatic (Spatially Correlated) 


Pulsar data are referred to an observatory time standard when 
they are produced, but they are transformed into Barycentric 
Dynamical Time (TDB) for timing calculations and compar- 
ison across data sets (Hobbs & Edwards 2012; Luo et al. 2021). 
Observatory clocks are not perfect, and so corrections are 
applied at various stages of this transformation. See Luo et al. 
(2021) for a detailed discussion of clock corrections and 
transformations used in pulsar timing. Any errors in these 
transformations, or drifts between time standards, can manifest 
in pulsar data as the same offset in all data sets with the clock 
error (Hobbs et al. 2012; Miles et al. 2023). These offsets are 
spatially correlated in the sense that all affected pulsars will be 
impacted by the same shift at the same time. The spatial 
correlation function is therefore monopolar, i.e., the same for 
all pulsar sky positions. See Tiburzi et al. (2016) for a complete 
discussion of clock errors and other spatially correlated noise 
processes in pulsar timing data. 


2.9. Measurements at the Observatories—Varying 
Chromaticity 


Radiometer noise from the observing systems used at our 
telescopes dominates our measurement uncertainty for our 
lowest-S/N pulsars but becomes a less significant component 
at higher S/N; we are radiometer noise dominated for most of 
our pulsars at most epochs. Random fluctuations due to 
radiometer noise should be time and frequency independent. 
However, the system temperature, and thus the S/N, is 
chromatic, as it depends on frequency-dependent contribu- 
tions from the receiver bandpass, Galactic background, etc. 
We have shown that the uncertainties from our template- 
matching procedure follow the expectations from radiometer 
noise for lower S/N (see Appendix B in Arzoumanian et al. 
2015 and note that we correct for this effect in our more 
recent data sets) and deviate owing to jitter and diffractive 
interstellar scintillation at the higher S/Ns (Lam et al. 2016), 
regardless of the frequency. Therefore, even though the 
system temperature will be frequency dependent, our per- 
channel modeling of the template-fitting error should contain 
no known biases. 

Incorrect polarization calibration could be a source of error 
in our data. For every observation we correct for differential 
gain and phase variations in the two chiralities of polarization 
of the received radio waves, as well as the changing parallactic 
angle. The latter effect is easily calculated based on the 
source's apparent position, whereas the former is corrected for 
based on measurements of a noise diode taken prior to each 
observation and an unpolarized quasar roughly once per month 
by which we correct for variations in the noise source power. 
Miscalibration causes alterations to the polarization-summed 
profiles used in timing (Stinebring et al. 1984), which would 
then be a chromatic source of noise given the frequency 
dependence of the (polarization) profiles. If polarization 
calibration is incomplete, for example, when a calibrator 
source is not observed at a particular epoch, we would expect 
some contribution to chromatic WN. 

Strong radio frequency interference (RFI) is removed in 
several stages of our pipeline as discussed in NG15. Low-level 
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RFI will still remain and perturb the TOA estimates in the 
template-matching procedure. RFI can take many forms— 
remaining narrowband RFI will perturb individual frequency 
channels, but broadband RFI can perturb all of the TOAs across 
the band. RFI can also be impulsive or periodic. Both can be 
removed from individual subintegrations if strong enough to 
identify. If not, in either case RFI will be folded at the pulse 
period and diminished, but it will still affect the TOAs as a 
source of possibly chromatic WN. 


3. Phenomenological Noise Model 


Independently fitting for each and every source of noise 
identified in Section 2 across all pulsars in the PTA would be 
extremely challenging. In order to efficiently characterize the 
noise of our detector and allow noise modeling from unknown 
sources, we therefore use a phenomenological model that takes 
these effects into account while reducing the number of fit 
parameters dramatically. The model has two main parts, 
distinguished by the timescale and spectral characteristics of 
the noise being modeled. The first part models the white, or 
uncorrelated in time, noise that has a "flat" contribution to the 
power spectra of timing residuals. We model this noise using 
three parameters that increase the uncertainties on the TOAs by 
accounting for WN unaccounted for in the template-fitting 
process (Arzoumanian et al. 2015). The second part uses 
Gaussian process regression (Williams & Rasmussen 2006) to 
model low-frequency, time-correlated, or red, stochastic 
processes in the data. In theory, both of these models could 
be included in the covariance matrix for the TOA data; 
however, it is more numerically expeditious to separate out the 
Gaussian process formalism for real data analyses. See 
Section 4 for details about how this is accomplished. 


3.1. White-noise Model 


The uncertainties that are initially associated with pulsar 
TOAs os,w are due to the finite S /N of the matched-filtering 
process used to calculate them (Taylor 1992). A template of the 
pulse shape, built from many observations, is convolved in the 
Fourier domain with the summed pulse profile from a single 
observing session, under the assumption that the data compose 
a scaled and shifted version of the template added to WN 
(Lommen & Demorest 2013). However, there are other sources 
of WN (e.g., jitter, scattering, RFI; see Section 2) that do not 
adhere to the assumptions of matched filtering (i.e., that the 
pulse shape is a copy of the template). These various noise 
effects cumulatively induce variance in the TOAs, and it is 
nontrivial to disentangle the distinct contributions of each noise 
source to the total TOA uncertainties. We therefore employ an 
empirical WN model for pulsar timing data that inflates the 
Os/n measured from the pulse template-matching process using 
three parameters. 

Three WN parameters are used to adjust the TOA 
uncertainties in order to accurately reflect WN present in the 
data. This process yields a reduced X? near unity for the fit to 
the timing residuals, if the timing model is complete and 
accurate. Various differences between pulsar timing back ends 
and radio observatory receivers make it necessary to give 
different values of these WN parameters to each receiver/back- 
end combination. These parameters thereby encode the 
trustworthiness of TOAs from each receiver/back-end combi- 
nation, down-weighting the TOAs from combinations where 
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effects in addition to template matching reduce the reliability of 
the data. These three WN terms—EFAC (F), EQUAD (Q), 
and ECORR (7)—come together with receiver/back-end 
combination re/be dependence as 


Cj = F(re/be) [og s; + Q°(re/be)] 6j + Jre[be) Ui, 
(2) 


where the i, j denote TOA indices across all observing epochs, 
bi is the Kronecker delta, and we omit the dependence on 
receiver and back end, re/be, from here on for simplicity. 
While EFAC and EQUAD only add to the diagonal of C, 
where Cj; are the elements of the covariance matrix to be 
discussed in Section 4, the ECORR terms are block diagonal 
for single observing epochs. ECORR is modeled using a block- 
diagonal matrix, U, with values of 1 for TOAs from the same 
observation and 0 for all other entries. 

Historically, the EFAC (error factor, F) parameter was the 
first WN parameter added to the pulsar timing covariance 
matrix. EFAC is a scale factor on the os/n. The increase in 
uncertainties from EFAC attempts to account for under- 
estimated template-matching errors from low-S/N observing 
epochs and template mismatches due to pulse profile 
variability. This parameter has tended toward 1 as pulsar back 
ends have improved, more high dynamic range (i.e., 8 or more 
bit sampling) digital systems have been implemented, and the 
treatment of low-S/N TOAs has improved (Alam et al. 
2021a, 2021b). 

To include additional WN, the EQUAD (Q) parameter?” 
(first used in Nice & Taylor 1995) is added in quadrature 
following the usual rules of noise propagation. EQUAD 
encompasses additional sources of WN not accounted for in 
the TOA uncertainties and not modeled by EFAC. 

The most recent WN parameter added to the pulsar 
covariance matrix is ECORR (error correlated in radio 
frequency, 7), which is also added in quadrature and 
specifically models noise correlated across radio frequencies 
within a single observing epoch.” ECORR was first used in 
Arzoumanian et al. (2015) to mitigate noise correlated across 
the narrowband TOAs obtained for a single observing epoch. 
In part this parameter became necessary owing to NANO- 
Grav’s data acquisition strategy of using multiple TOAs across 
the full observing band for chromatic mitigation. While 
ECORR is not strictly “WN” from the perspective of being 
only diagonal in the covariance matrix, its spectral contribution 
is effectively white, i.e., constant in the frequency domain, for 
the frequencies important to GW searches. The correlation 
timescale modeled is only across pulses emitted from the pulsar 
within a relatively short (~30 minutes) observation. ECORR 
partially accounts for intraband correlations caused by pulse 
jitter, but it also includes other short-timescale correlations 
across radio frequency, which can be produced by, e.g., RFI 
and time-variable scattering (Lam et al. 2019; Shapiro-Albert 
et al. 2021). 


92 Note that there are two distinct definitions of EQUAD in the literature, 
depending on whether EFAC multiplies only og/n, referred to as the 
TEMPONEST convention for the paper (Lentati et al. 2014) where it was first 
used, or whether it multiplies the sum, in quadrature of osx and EQUAD. All 
three of the main pulsar timing packages and the Enterprise software stack 
use the latter convention, laid out in Equation (2). 

70 Here an observing epoch is effectively one observation with a single 
receiver. Multiple observations within a single MJD are not treated as 
correlated. 


THE ASTROPHYSICAL JOURNAL LETTERS, 951:L10 (57pp), 2023 July 1 


3.2. Red-noise Model 


Various sources of noise in Section 2, such as variations in 
the electron density of the ISM, instabilities in the spin of the 
pulsars, corrections to Earth-based clock systems, corrections 
to the solar system ephemeris model, and, of course, TOA 
shifts due to the stochastic GWB are time correlated across 
long timescales. The characterization of these effects in the 
frequency domain, e.g., a power spectral density, shows that 
the effect has more power at lower (redder) frequencies, hence 
RN. In most cases the theoretical models for these sources of 
noise have a power spectral density P(f) that follows a power 
law with frequency f (Cordes & Shannon 2010), P(f) ~ A?f ^, 
with amplitude A and where the “redness” is dictated by the 
explicit minus sign in the exponent and ^ > 0. More complex 
models have been investigated, such as turnover models, but 
the straight power law is still considered an accurate and 
effective model for MSPs (Goncharov et al. 2020). 

There are two methods by which RN can be included in a 
pulsar timing analysis. One way is to use the Wiener—Khinchin 
theorem to write the correlations between different TOAs, 
(ti, tj). 


T 
=f df cos[2nf (t; — t] PC). (3) 


where CAN are the elements of the covariance matrix as noted 
above and fwy4 is the Nyquist frequency. Note that, strictly 
speaking, the integral in Equation (3) does not converge for all 
forms of P(f), in which case a low-frequency cutoff must be 
used; see van Haasteren et al. (2009) for an exhaustive 
discussion. The details of the RN model are discussed more in 
Section 4. 

Alternatively, the perturbations of the TOAs due to RN can 
be modeled directly with a Fourier basis and a set of 
coefficients, Fa, where F is an (Nroa X 2Nfreq) Fourier design 
matrix and a is a 2Ngeg vector with two coefficients for each of 
Nereq frequencies included. As will be discussed in Section 4, 
the form of P(f) can be dictated by imposing various functional 
forms on the a coefficients. This method is at the heart of PTA 
searches for the GWB (Lentati et al. 2013; van Haasteren & 
Vallisneri 2014; Lentati et al. 2016) and achromatic (in radio 
frequency) RN models used for individual pulsars. In practice, 
the methods of Gaussian process regression are used in the 
analyses to accurately include the stochastic nature of these 
signals, and the explicit Fourier basis modeling will eventually 
be incorporated into the covariance matrix, but in a different 
form; see Section 4. It is these analyses that are used to fit for 
RN parameters used in the covariance matrix for pulsar timing 
packages. While NANOGrav uses the DMX model (NG15) in 
this analysis, these same Gaussian methods are used by other 
PTA collaborations to model the variations in DM by including 
a dependence on the radio frequency (Lentati et al. 2016). 


4. The NANOGrav Pulsar Covariance Matrix 


PTA data analysis is done using the timing residuals ôt, 
produced by subtracting TOAs predicted by a timing model ty 
from the TOAs t: 


ôt = t — ty = tp — ty + Fa +n, (4) 


where the boldface symbols are column vectors. The difference 
between the actual underlying, deterministic, non-GW delays 
and the timing model, tp — tmy = Me, is represented by a linear- 
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order Taylor series in the (assumed small) perturbations to the 
model parameters, óp — e, and the design matrix, M, where 


Ot i 
m= (2 (5) 
Pj P=Po 


is evaluated (usually analytically, numerically otherwise) at the 
best-fit parameters, po, from a linear least-squares analysis. The 
use of the linearized timing models was introduced in Ellis 
et al. (2013) and is especially expeditious for full PTA 
analyses. However, full timing model fits can also be done as 
part of the Bayesian searches (Lentati et al. 2014; Vigeland & 
Vallisneri 2014, A. R. Kaiser et al. 2023, in preparation). The 
Fa term encodes the long-timescale correlated (red) noise, for 
both individual pulsar data sets and the GWB; see Section 3. 
Finally, n is the WN remaining in the residuals, assumed to 
have a multivariate Gaussian distribution, 


exp(—4n"C-'n) 
Jdet2nC) 


with covariance C(E) = (n, n”), which is a function of the set 
of WN parameters, e.g., E = (F, Q, J). We construct the 
likelihood function for a PTA by combining the timing 
residuals, C, and Equation (6): 


P(n) = (6) 


exp( 6t — Me — Fa) C Wét — Me — Fa)) 


Jdet(27C) 


P(étle, a, E) = 


(7) 


When conducting a search for GWs in pulsar timing data, we 
are not concerned with the values of the timing parameter 
perturbations included in e. Using a Gaussian distribution for 
the prior probability distribution of the parameter perturbations, 


exp(—ye"Xx-le) 
Jdet(27X) 


and defining X = (e, E. one can marginalize over the timing 
model parameter perturbations (see Appendix A.1) 


P(e) = (8) 


P(6t\a, E) = f T(6t|e, a, E)P(e)d(e) (9) 


exp( — ; 8t — Fa)’ D- (êt — Fa)) 
Oe (10) 
det (27D) 
where D = C + MXM” (Equation (A3)). 
In a similar fashion the Fourier basis coefficients are also 
marginalized over assuming that the RN, whether intrinsic to 
the pulsar or a GWB, is a Gaussian process, 


exp(—4a" pa) 
/det(27y) 


where y= (a, a’) is a (2Nfreq X 2Nfreq) diagonal matrix with 
dimensions dictated by the number of frequencies, Ngeg, used 
in the analysis and describing the variance of the RN 
coefficients. In principle, these coefficients could be used as 
free parameters; however, they are usually parameterized as a 


Pa) = 


(11) 
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Table 1 
Connecting Noise Terms: Astrophysical Noise Sources, Along with a Summary of Their Attributes and Connections to the Phenomenological Noise Model 

Noise Type Symbol Origin Phenomenological Model Component Timescale Spatially Correlated 
Radiometer Os/N Telescope OS/N Short 

Jitter 0j Pulsar I Short 

Diffractive interstellar scintillation OpIss ISM Q, JT Short 

DM misestimation óDM ISM RN/DMX Long 

Solar wind DM misestimation óDM,, IPM RN/DMX Long F4 
Frequency-dependent DM DM(v) ISM RN/DMX Long 

RFI Telescope I Short 

Polarization miscalibration Telescope T Both 

Scattering ISM dU Short 

Solar system barycenter mismodeling solar system RN Long "4 
Clock errors Telescope RN Long "4 


Note. Short timescales are roughly the order of the observation length or less, while long timescales refer to periods between observation epochs. 


power law, account for the noise processes of the NANOGrav-observed 
pulsars; see Table 1. Collectively, the WN in each of our pulsar 
zs -7 Ag f =y timing models can be expressed as a covariance matrix where 
inis ; Es yr?df ~ ——L—— 3t yr, (12) each of Nroa TOAs is specified with both its recorded time t 
127° 1 fs, 1207 Ta \ yr and its radio frequency v. This covariance matrix is then fit 
globally over all epochs and over all of a pulsar's front-end/ 
or some other functional form for the power spectral density. back-end combinations. For a set of TOAs on a given epoch, 
Note that fj, — 1 yr !, and the “~” assumes the usual constant we write the noise model in the form 
Flog nv) + Qe + 7? I? e J? 
epoch — g’ Flag NU) + Q? - J? -- J? | (14) 
g? I? + Fossa) + Q9 + I? 


df~1/ Tspan, Where Tspan is the total time span (years) of the where F, Q, and J are the EFAC, EQUAD, and ECORR 


data set. Marginalizing over the a coefficients gives the final described in Section 3, respectively. This matrix is block 
form for a single-pulsar likelihood, diagonal, with the entries of J? connecting frequency channels 
from the same observations in blocks. The WN covariance 
PAE) = [re étla, E)P(a)d(a) matrix can then be written down as 
epoch "T 

exp( — 16t7N- ét) Cop) j o 

= ; (13) E 0 CP (h) + 0 
Jdet(2TN) C= : 7 (15) 

epoch 

where N= D + FoF”, again see Appendix A.1 (Equation (A7)). 0 0 se CPP) 
In a single-pulsar noise analysis using this likelihood, one 

searches over all of the parameters (A, y, E); however, the WN where the noise is independent between epochs, but the noise 
parameters are fixed to their maximum likelihood values parameters are fit globally over all TOAs from a single pulsar. 
during full PTA GW searches. The WN parameters are If one wants to include the RN in the covariance matrix, the 
independent of the parameters varied in the GW search, and total noise would then be C — C+ CP". This approach is used 


the necessary matrix inversion’! is prohibitively expensive 
from a computational standpoint. The marginalizations above 
were developed to reduce this computational expense to a 
single inversion of a (2M freq x 2Nfreq) Matrix (van Haasteren & 
Vallisneri 2014). 


in pulsar timing software when the best-fit amplitude and 
spectral index for RN in a given pulsar have already been fit for 
in an earlier analysis (Coles et al. 2011). 

It is tempting to try to map the phenomenological noise 
model to a signal covariance matrix made up of all of the 
effects from Section 2. However, in practice, there are a few 


Sel EAA RENOIR CONG EE MAGI reasons for not doing so in a GW analysis. First of all, the 


As noted above, the phenomenological WN model is a complexity of these models makes this task unwieldy for GW 
function of the parameters E = (F, Q, J), which together searches, especially given that some of these effects are not 
directly measurable. Second, the noise measured through the 

7! For further mathematical discussion of this inversion, see Appendix A.1. phenomenological model is known to be larger than the sum of 
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the various physical effects, so a signal-based covariance 
matrix might not include all of the noise in our detector. In 
Section A.2 we sketch a signal covariance matrix to show how 
some of the terms from Section 2 would enter into the analysis. 


5. PTA Noise Budget 
5.1. Comparison of Physical and Phenomenological Models 


The noise in Section 2 enters into the phenomenological 
noise model through a diverse set of physical processes, and 
any particular physical effect may enter more than one part of 
the phenomenological model. From the standpoint of detecting 
GWs with PTAs, the most salient aspect of this complex 
relationship is that the phenomenological model includes the 
known physical sources of noise and is pliable enough to 
include unknown sources as well. 

That pliability is demonstrated in Table 1, where we provide 
a list of the various noise sources from Section 2 and 
relationships between the phenomenological parameters in 
which they are included, along with a number of other 
descriptors that can be used to sort the sources of noise. 


5.1.1. DM Misestimation due to Asynchronous Measurements 


In order to estimate DM on each epoch, we make TOA 
measurements across a wide frequency range. We used two or 
more receiver bands at one or more telescopes to cover the total 
range. Data were recorded with each receiver independently, 
with observations occurring within 30 minutes of each other at 
Arecibo but within several days at Green Bank owing to 
efficiency constraints in changing the receivers. DMs can 
change rapidly (i.e., on timescale of a day or less) owing to the 
passage of the line of sight through the inhomogeneous ISM, 
changes in the solar wind, etc. In the presence of this stochastic 
variation, Lam et al. (2015) estimated the effect of asynchro- 
nous measurement, finding an induced perturbation on the 
infinite-frequency arrival time with a shallow (y= 2/3 for a 
Kolmogorov medium) RN spectral index. The rms error is 
given as (see Equation (14) in Lam et al. 2015) 


t, sep 


v \? r? 3/6 
Og... Z 6.5 nd ) . 
GHz r? — 1A Ata: Guz/1000 s 


(16) 


This expression assumes a simplified model where pulsar 
observations are made at only two spot frequencies v and v, 
where r = v/v' with v > v’, tsep is the time spacing between 
the observations, and At, 1 Guz is the scintillation timescale at a 
radio frequency of 1 GHz. For pulsars observed at Green Bank 
with r= 1.5 GHz/0.8 GHz = 1.875, a time between observa- 
tions (in days) of fseẹp 7 2, and a typical scintillation timescale 
of approximately 1000s at 1 GHz, we see that the rms error 
from this effect is approximately 7 ns. These errors are 
therefore generally smaller than the intrinsic RN we measure 
for most pulsars. For PSR J17134-0747, observed with all three 
telescopes, reducing the impact of rapidly varying DM requires 
the window over which we fit a single DM to be shorter than 
typical (e.g., <1 day compared to 6 days) to avoid a larger DM 
misestimation error. 
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5.1.2. DM Misestimation due to Additional Chromatic Effects 


There are a host of unmitigated chromatic effects in our 
timing model, though many or all are expected to be small 
except for specific lines of sight. Scattering is expected to be 
the next dominant effect. It is highly frequency dependent, 
broadening the width of the PBF of the ISM by an amount that 
is xv ^^ for a Kolmogorov medium. When this scattering 
delay tscat K W, where W is the pulse width, the first-order 
effect is to delay the TOA by tcan. However, when tscatt is 
comparable to W, there is distortion of the pulse shape and the 
time delay becomes a nonlinear function of the scattering delay 
(Hemberger & Stinebring 2008). We are likely only in this 
regime for one NANOGrav pulsar, J1903+-0327, with a spin 
period of 2.2 ms and scattering delay of ~150 us at 1.5 GHz 
(Geiger & Lam 2022). 

As discussed in Section 2.5, DM itself is a frequency- 
dependent quantity (Cordes et al. 2016). Again considering a 
simplified model where pulsar observations are made at only 
two spot frequencies v and v’, we can write the TOA 
perturbation as 


Fico 


Icy! + Ey! 


Stoo = ; (17) 
where fc, and tc," are nondispersive chromatic delays at 
frequencies v and v’, respectively, and e, and c, are additive, 
frequency-dependent errors at frequencies v and v’, respec- 
tively. The additional error is uncorrelated between different 
frequencies for radiometer noise but is highly correlated for 
intrinsic pulse jitter, at least over modest frequency separations 
(see below). The rms timing error is found by calculating the 
square root of the variance. Even when there is no WN, i.e., 
€, = & = 0, there is a nonzero perturbation that affects our 
timing. Frequency-dependent DM affects us at the tens of 
nanosecond level (Lam et al. 2018b) for most pulsars at the 
frequencies we observe, though for a few at the highest DMs 
that number can be substantially higher. 


5.1.3. Accounting for Jitter 


As noted in Section 3, the ECORR parameter was added to 
the pulsar signal covariance matrix to account for frequency- 
correlated WN. When examined over NANOGrav’s 12.5 yr 
data set pulsars, we find that a global fiducial estimate for the 
rms uncertainty from jitter is (Lam et al. 2019) 


60 w( B | Tops ) 
Oj ~ 
: 5 ms 30 minutes 


for a pulsar with a P, — 5 ms spin period and an observing time 
of Ty, = 30 minutes, if we take the 1500 MHz estimate as a 
representative value. This uncertainty level makes up a 
nonnegligible fraction of our known WN budget. In its general 
use, for narrowband timing, the ECORR estimation is higher 
than the physically derived constraints in Lam et al. (2019; see 
also the comparison between ECORR and jitter in Lam et al. 
2016), which prevents potential underestimation of the 
frequency-correlated noise. Likewise, the in-depth study of 
jitter in PTA pulsars revealed that jitter typically decorrelates 
given a sufficiently broad frequency range (Shannon et al. 
2014). This decorrelation bandwidth is currently not incorpo- 
rated into the ECORR parameter because current receivers are 
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narrow enough to avoid decorrelation; however, with the 
implementation of wideband receivers in the near future (see 
Section 6), this decorrelation bandwidth may fall within the 
band of these new receivers. 


5.1.4. The Impact of Polarization Miscalibration 


For any given narrowband profile, polarization gain 
miscalibration will cause a TOA uncertainty dependent on 
the fractional gain error € and the degree of circular 
polarization 7y for a pulsar with pulse width W as (Cordes 
et al. 2004) 


bords (E) Ww 
Las TP CEP LETTRE I 


The fiducial values given are representative of the NANOGrav 
pulsars and data (Lam et al. 2018b), though the my and W 
values can still vary quite widely (e.g., van Straten 2006; 
Gentile et al. 2018). Cross-coupling induces a false circular 
polarization that adds to the measured circular polarization. 
While a seemingly large value, the uncertainty is on the 
narrowband TOAs and not the infinite-frequency arrival time 
(e.g., the epoch-averaged TOAs), and so the value should be 
reduced by the square root of the number of TOAs, a reduction 
of nearly an order of magnitude. Polarization miscalibration 
uncertainties in the NANOGrav data set have not been fully 
explored but are expected to add of order 100 ns uncertainty to 
our excess noise budget. 


(19) 


5.2. Red-noise Analysis 


Since the strongest GW signal expected in PTAs manifests 
as a long-time-correlated stochastic signal, understanding RN 
in pulsar timing data sets is of the utmost importance for 
mitigating false positives and/or parameter misestimation. The 
use of power-law RN models in pulsar timing analyses is 
ubiquitous and predates modern PTA data analysis techniques 
(Blandford et al. 1984; Cordes & Downs 1985). There are a 
number of more complex models for RN in PTA data sets (Lam 
et al. 2018a; Hazboun et al. 2020b; Goncharov et al. 2021a; 
Chalumeau et al. 2022), but we will demonstrate shortly that 
the power-law model is a reasonable and effective choice by 
carrying out an additional noise analysis. 


5.2.1. Accounting for Red Noise 


Searches for GWs in PTA data sets start with single-pulsar 
noise analyses and then move on to full analyses that include 
all of the pulsar data sets. In the case of the GWB this involves 
a correlated search across all pulsar data sets. The main reason 
for the individual noise analyses is to search for the WN 
parameters described in Section 3. They involve full 
Nroa X NqoA (O(10*) narrowband TOAs per pulsar, and a 
factor of ~30 reduction on average for the wideband data) 
matrix inversions and upward of 21 parameters per pulsar. 
Hence, they are too costly to do across all pulsars at once in a 
full Bayesian search. It was shown in early, smaller data sets 
(Demorest et al. 2013; Arzoumanian et al. 2015) that the WN 
parameters do not change substantially when the analysis is 
done across other pulsars, as long as there is an RN model 
included in both analyses, since this model can substantially 
change the amount of WN in some pulsars. The presence of 
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Table 2 
Individual RN Model Parameter Values and 68% Credible Intervals for Pulsars 
with Significant Detections of RN 


Pulsar log) ARN YRN 

B18554.09 — 14.0423 3.9708 
B1937421 —13.6+2! 4.0703 
B19534-29 12.802 18°93 
J0030+0451 —14.4+04 4.6109 
J0437—4715 — 13.4403 0.5706 
J0610—2100 —12.9*03 44*22 
J0613—0200 —13.8/03 3.107 
J1012+5307 12.610 0.8703 
J1600—3053 —13.5792 1.6553 
J1614—2230 —14.9*12 4.7735 
J1643—1224 -12351 0.90 
J1705—1903 -12.6701 0.503 
J171340747 —14.101 2.603 
173810333 —14.6*08 52413 
J1744—1134 —14.1*04 3.6*15 
J17454-1017 —119701 24708 
J1747—4036 -12.6*04 24449 
J1802—2124 —12.2403 18:04 
J1853+1303 -13.5:92 2.3495 
3190340327 -122:0! 1501 
J1909—3744 -145'01 4.1795 
J1918—0642 —13.gt04 2.7419 
J1946+3417 -12.550] 140 
J2145—0750 —12,940! 0.6133 
J223440611 —13.9*03 32413 


Note. Pulsars with italicized entries only have significant detections in 
individual pulsar noise analyses, while boldface entries show RN even when a 
common process is included. Note that Agw is in strain amplitude to match the 
GW analyses; hence, it is unitless. 


two RN analyses, one for individual pulsars and one in the full 
GWB search, allows us to do an accounting of the RN in each 
pulsar and keep track of which pulsars the RN power moves 
from the individual noise channel into the common signal. 

In Table 2 we show the pulsars that have significant RN 
detections in their individual noise analyses. Our criterion is the 
same as for NG15; we consider the detection to be significant 
when the pulsar's RN has a Bayes factor 6 > 100, using the 
Savage—Dickey (Dickey 1971) approximation when possible. 
The italicized entries in Table 2 denote the pulsars that no 
longer have significant detections of RN in the full PTA 
analysis, i.e., the power in the single-pulsar RN model has 
moved into the common channel." The boldface entries denote 
the pulsars that continue to have significant RN detections. 

Figure 2 shows the RN power-law posteriors for the same 
pulsars included in Table 2. Those italicized in Table 2 are 
shown in gray in Figure 2, along with the overlaid posterior for 
a common process across all of the pulsars. Note that the gray 
contours follow a well-known covariance trend between the 
spectral index and RN amplitude. See NGI5gwb for more 


72 Depending on the search, this can be either a spatially correlated process, as 
in searches for the GWB, or just a “common” process, described by only the 
power spectral density and not the spatial correlations. See NG15gwb for more 
details. 
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Figure 2. Two-dimensional parameter posterior contours for the power-law RN 
models used in the pulsar noise models of the NGI15 data set. The power 
spectral density is parameterized by a spectral index and amplitude. The black 
contours show the RN parameters of pulsars that have significant individual 
RN in both their single-pulsar analyses and the full PTA analysis. The gray 
contours show RN parameters of pulsars that have significant individual RN 
detections only in the single-pulsar analyses, i.e., the power in the individual 
pulsar model has moved into the common model. The blue concentric contours 
show the parameter recovery from a full PTA common process analysis 
(NGI5gwb). Posteriors for pulsars without significant individual RN detections 
are not included for clarity. 


details, along with a way to ameliorate the covariance in 
parameter recoveries. 

The results are broadly as one would expect. The pulsars that 
possess significant RN detections in their individual noise 
analysis, which lie along a similar slope of covariance to the 
common process recovery, lose their RN power to the common 
process. This is corroborated by a "leave one out" analysis 
in NGI5gwb, where individual pulsars are dropped out of a 
common uncorrelated red process (CURN) analysis and an odds 
ratio is calculated between an analysis with N pulsars and one 
with N — 1 pulsars. All of the grayed-out pulsars have positive 
dropout factors, i.e., the analysis with that pulsar is favored 
over the analysis without that pulsar, except for PSR J0437 
—4715 and PSR J1713--0747. These two pulsars are known to 
have challenging noise properties that are discussed below. 

It is important to note that while there is obviously another 
intriguing cluster of posteriors that are a few orders of 
magnitude larger than the recovered CURN, the RN in those 
pulsars does not seem to be common and would be excluded by 
the pulsars in the gray cluster that do not see this RN. Various 
Bayesian analyses find no CURN in that part of parameter space 
(NG15gwb). Given the shallow spectral index, this RN is likely 
due to mismodeled chromatic noise (Cordes & Shannon 2010). 
While only the individual noise analysis results are shown in 
Table 2, there is no change in the median RN values at the 
precision reported and only slight changes in the credible 
intervals. 

The individual RN models are working as expected. For 
pulsars where the RN seems to be a part of the CURN the RN 
models are a stand-in for the common process in the individual 
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noise analyses but do not hold onto the RN during the full PTA 
analysis, where the common process instead models the RN. 
For pulsars where the RN seems to be truly intrinsic to only a 
single-pulsar data set, the RN models are successfully modeling 
that noise, which is not able to contaminate the common signal. 
One interesting characteristic of the grayed-out posteriors is 
that they cover a large span in the spectral index, but all follow 
a roughly common slope. Pulsars with short time spans, e.g., 
PSR J0437—4715 with only 4.8yr of data and the gray 
posterior against the left side of Figure 2, may have shallow 
spectral index recoveries because they lack enough sensitivity 
at low frequencies where the GWB will rise above their 
WN floor. 

Next, we present three separate noise analyses for four 
separate representative pulsars^ in Figures 3—6. Each figure 
contains all of the components of the power-law model, a free 
spectral analysis, and what we refer to as an "excess noise 
analysis," discussed below. The plots show rms fluctuations of 
the timing residuals, i.e., the amplitude of the noise in the 
residuals across the frequencies shown, in units of log, (5). We 
convert the parameter posterior probabilities from all noise 
models into amplitude for ease of comparison against the 
precision of the measurements. The most generic noise model 
uses a separate set of y coefficients (Equation (12)), allowed to 
vary freely, without reference to any particular power spectral 
density model. This so-called free spectral model (Lentati et al. 
2016; Hazboun et al. 2020b) is a Bayesian spectrogram of the 
pulsar data. Here the spectra have been recovered from each 
pulsar's data across 30 frequencies, ranging from 1/T;,4, to 
30/ Tepan, Where Tspan is the full time span of the PTA.” The 
solid orange “violin” plots show the posteriors, frequency by 
frequency, for the free spectral analysis. Significant posteriors 
are those where the violin plot is separated completely from the 
frequency axis, as in the second-to-lowest-frequency free 
spectral (solid orange) posterior of Figure 3. Insignificant 
posteriors have broad posteriors all the way down to the 
minimum values, while slightly significant detections have thin 
posteriors extending to the minimum. 

The free spectral model posteriors represent a “raw 
spectrum, which does not assume anything about the power in 
a pulsar data set However, when we search for common 
uncorrelated processes across the pulsars, a simpler power-law 
model is used for both the RN and the CURN. The recovered 
CURN from NG15gwb, plotted the same for all pulsars, is shown 
as the dashed red line. The horizontal gray dashed line shows the 
WN power spectral density of the residuals (subscript R), 


converted into these units as WN = ESI Tipan- The 
individual RN is shown with a dotted green line. This line 


may not be visible for pulsars where the power-law amplitude is 
small or not very significant. Since an individual RN model is 
used for every pulsar, regardless of significance, we have 
included the maximum likelihood values of the RN for all 
pulsars. We do not include the RN from individual pulsar noise 


» 


73 The choice of priors is very important to ensure that this is the case 
(Hazboun et al. 2020a). 

Ja Figures and analyses are included for all pulsars at the end of this 
manuscript. 


75 The choice of these 30 frequencies is based on the usual Nyquist frequency 
sampling considerations and extends up to frequencies higher than where a 
GWB would be detectable owing to the WN floor (Lentati et al. 2013). 
Unfortunately, the use of these models, with 30 frequency parameters for each 
pulsar, is currently not feasible for full PTA analyses owing to computational 
limitations. 
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NG15yr Noise Budget: J1909—3744 (Typan=15.5 yrs) 
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Figure 3. The excess timing residual delay as a function of frequency for PSR J1909—3744. The total noise (solid black line) includes WN (gray dashed line), 
common RN (red dashed line), and individual RN (green dashed line). The free spectral model (solid orange posteriors) does not dictate any relationship between the 
amplitude of the power at different frequencies. Holding the total noise model (solid black line) parameters constant and searching for additional noise using a free 
spectral model results in the excess noise shown in blue. The vertical dashed line denotes a frequency of 1 yr |. Bayes factors for these parameters, shown in the 


bottom panel, are fairly insignificant across all frequencies. 


analyses in these figures for clarity. See Table 2 below for a list 
of pulsars with significant detections of RN. The solid black line 
shows the total of the common process RN, the WN, and any 
individual RN model and represents the total noise power (in 
these units) for the pulsar using these models. 

Sources of RN were studied in Cordes & Shannon (2010). 
Large RN with a shallower spectral index is thought to be due 
to modeling errors of time-correlated chromatic effects, while 
steeper spectral indices are thought to originate from achro- 
matic processes either intrinsic to the pulsar system or due to a 
stochastic GWB. The pulsars can broadly be separated into four 
categories depending on what type of noise they are dominated 
by in their lowest frequencies: 


1. Common Process Dominated: Figure 3 shows the noise 
budget for PSR J1909—3744, one of the most sensitive 
pulsars in the NANOGrav array. This pulsar has very low 
WN power spectral density, has a significant detection of 
power-law RN in its individual noise analysis, and shows 
significant recovery of power at a few low frequencies of 
the free spectral noise analysis. As denoted by the 
italicized text in Table 2, most of the RN power moves 
into the CURN in a full PTA analysis; hence, it is 
dominated by RN that appears to be a part of the common 
process. See Section 5.2.1 for a full accounting of the RN 
across the PTA. Looking at the free spectral posteriors, 
one might ask whether a power-law model is sufficient for 
such a pulsar, but in fact the excess noise analysis shows 
no significant detections of additional noise. The standard 
model accounts for power in the 30 frequencies 
considered here. This is true across all of the pulsars in 
the array. 
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2. WN Dominated: Figure 4 shows a representative of a 


second category of pulsar—one dominated by WN at the 
lowest frequencies. PSR J0509--0856 only has 3.6 yr of 
data, just making the cut for a pulsar included in 
NANOGrav gravitational analyses. One does not expect 
sensitivity to RN from this pulsar at the lowest frequencies 
considered here, since they are based on the longest time 
span covered by the full PTA data set. One can see that 
there are no significant detections of power at individual 
frequencies in either the free spectral model or the excess 
noise model. 


. Shallow RN Dominated: Figure 5 shows a third type of 


pulsar, PSR J19034-0327, also dominated by RN at the 
lowest frequencies, but which does not move into the 
common channel during a full PTA analysis. PSR J1903 
+0327 is our highest DM pulsar, and so significant 
unmodeled chromatic propagation effects are expected to 
impact the timing (see, e.g., Geiger & Lam 2022). Here 
again, the power-law RN model is effective at modeling 
most of the noise picked up by the free spectral model, as 
seen by the insignificant Bayes factors. 


. Steep RN Dominated: Figure 6 shows a fourth type of 


pulsar, also dominated by RN at the lowest frequencies, 
that also does not move into the common channel during a 
full PTA analysis. PSR B19374-21 is a well-known 
example of high-amplitude, steep spectral index RN that 
appears to be intrinsic to the pulsar system. Here again, 
the power-law RN model is effective at modeling most of 
the noise picked up by the free spectral model, as seen by 
the insignificant Bayes factors. 
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NG15yr Noise Budget: J05094-0856 (1.544, —3.6 yrs) 
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Figure 4. The excess timing residual delay as a function of frequency for PSR J0509+-0856. See Figure 3 for details. Bayes factors for these parameters, shown in the 
bottom panel, are insignificant across all frequencies. WN dominates this pulsar's short-time-span data set. 


NG15yr Noise Budget: J19034-0327 (1.54, —10.7 yrs) 
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Figure 5. The excess timing residual delay as a function of frequency for PSR J1903+0327. See Figure 3 for details. Bayes factors for these parameters, shown in the 
bottom panel, are fairly insignificant across all frequencies. The lowest frequencies are dominated by RN intrinsic to the pulsar. 


Noise budget figures for all of the pulsars used in Agazie run an additional noise analysis on each pulsar. This analysis 
et al. (2023b) are included as Figures 11—73. fixes the WN, power-law RN, and GWB parameters, i.e., the 
. solid black line in Figures 3 and 4, and adds in an additional 

5.2.2. Comparison of Time-correlated Noise Models analysis, using a separate set of y coefficients, allowed to vary 

In order to understand the quantities of noise that might be freely. This additional model is effectively the same as the free 
unaccounted for by the power-law model in our data set, we spectral analysis mentioned above, but this model now 
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B1937+21 (Tspan=15.9 yrs) 
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Figure 6. The excess timing residual delay as a function of frequency for PSR B1937+21. See Figure 3 for details. Bayes factors for these parameters, shown in the 


bottom panel, are fairly insignificant across all frequencies. The lowest frequencies 
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Figure 7. Free spectral (orange) and excess noise (blue) model Bayes factors for all 67 pulsars used in the 15 yr GWB search. The orange circles with arrows represent 
posteriors where the detection is too significant to use the Savage—Dickey approximation. The blue and orange circles represent the same frequencies, but the blue 
circles are offset slightly to the left of the orange for ease of viewing. The vertical dashed black line demarcates the 14th frequency. Only the lowest 14 frequencies 


were used in the GWB search. 


accounts for any excess noise unaccounted for in the usual 
noise model. As we will see, it is mostly superfluous when the 
power-law model is used. The WN parameters are from the 
individual pulsar noise analyses discussed in Section 3.1. In 
addition, we use the two-dimensional posterior power-law RN 
maximum likelihood values from a full PTA search for the 
CURN, along with the median value for the CURN amplitude from 
that same analysis. The objective is to understand the noise 
unmodeled by the WN, RN, and CURN stochastic processes. As 
with the free spectral model, the analyses are done over 30 
frequencies ranging from 1/ Topan to 30/ Tspan The new free 
spectral coefficients, referred to as “excess noise,” quantify how 
much power at various frequencies in each pulsar’s data set is not 
modeled by the WN + RN + CURN model. 

The results from this excess noise analysis are compared to 
the results of the other noise analyses for the four representative 
pulsars in Figures 3—6. The excess noise posteriors are shown 
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as the open blue violins. The bottom panel shows the Bayes 
factor, calculated using a Savage—Dickey approximation 
(Dickey 1971), for the excess noise parameters. 

Figure 7 summarizes the excess noise results from all of the 
67 pulsars in the NANOGrav PTA. We show the individual 
Bayes factors for each pulsar across the 30 frequencies 
considered in our GWB search. Again, the filled orange circles 
show the Bayes factor for the free spectral parameters, while 
the blue circles show the Bayes factors for the excess noise 
parameters. The orange circles with arrows at the top represent 
posteriors where the detection is too significant to use the 
Savage—Dickey approximation. The number of such detections 
is noted with a digit in the circle. The main takeaway from 
Figure 7 is that there is no case where a significant free spectral 
detection corresponds to a significant excess noise detection, 
revealing that the power-law model used across the PTA is a 
sufficient model, given the sensitivity of the data, for mitigating 
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RN in individual pulsars and detecting RN from a common 
process. 


5.3. Sensitivity Curves 


Detection sensitivity curves are commonly used in the GW 
community to summarize all aspects of detector characteriza- 
tion into a single "figure of merit." They are often used by the 
broader astrophysics community to assess the detectability of 
various GW sources. Here we use the Python package 
hasasia (Hazboun et al. 2019a), based on the formalism 
developed in Hazboun et al. (2019b), to calculate sensitivity 
curves for the 67 pulsars in the NANOGrav PTA and then 
combine these into a sensitivity curve for the stochastic GWB. 
This combination is based on the S/N of the GWB optimal 
statistic developed in Anholm et al. (2009), Chamberlin et al. 
(2015), and Rosado et al. (2015). 

The effects of the timing model on searches for various 
astrophysically interesting signals, especially noise correlated 
on long timescales, have been known for a long time 
(Blandford et al. 1984). These effects can be represented by 
a transmission function, 7(f), that encodes the fraction of 
power transmitted through the timing model. In this way the 
transmission function acts as the transfer function for a pulsar 
(Blandford et al. 1984; Cordes & Downs 1985). A slightly 
different, but equivalent, method (van Haasteren & 
Levin 2013) than described in Section 4 is used to margin- 
alize over the timing model parameters and build the 
transmission function. Here the G matrix, derived from the 
timing model design matrix, M, encodes the information 
about the timing model fit and acts to project the data and 
covariance matrix into a basis orthogonal to the timing model 
(see Hazboun et al. 2019b for details). The transmission 
function can be written in terms of G, 


1 NO (GGT )je??nf in, 
TOA ij 


KA= (20) 


and calculated for any pulsar in the array. In Figure 8(a) the 
transmission function for PSRJ1909—3744 is shown. The 
transmission function has a few interesting features. First, the 
fit for the rotation frequency and its derivative (the spin and 
spin-down) of the pulsar pulls a quadratic polynomial out of the 
data, which acts as a high pass filter, and limits sensitivity at the 
lowest frequencies. As the time span of the pulsar data set 
increases, the spin-down parameters are fit more precisely and 
the frequency at which the transmission “turns down" moves to 
lower and lower frequencies. Second, the single-frequency fits 
for the sky position/proper motion of the pulsar (f= lyr”) 
and parallax (f— 2 yr ') remove power in a narrow band 
around those frequencies. If the pulsar has a binary period 
within the frequency range, there will be another dip in the 
transmission function for that fit as well. The width of these 
dips is proportional to 1/Tspan and therefore narrows the 
longer a pulsar is timed. Lastly, the DM variation model used 
removes power across the entire GW frequency band. The 
DMX model constitutes >200 parameters for a few of the 
pulsars timed and diminishes the power by a factor of ~3 
across the frequencies searched for GWs. Therefore, the 
transmission function does not asymptote to 7( f) ~ 1 at the 
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Figure 8. Various spectra for PSR J1909—3744. The transmission function is 
shown in panel (a) as a dashed orange line, while perfect transmission 
(Kf) = 1) is shown as the dotted black line. The various frequency-dependent 
features are discussed in the text. Various power spectral densities of the 
residuals, P&( f), are shown in panel (b). The black dotted line shows the WN 
power spectral density PON), The orange dashed line shows PY with the 
effects of the timing model (TM) included. The gray dashed line shows the 
GWB reproduced using the median values of the CURN” model recovered 
in NG15gwb. The blue solid line shows the full power spectral density of the 
pulsar and includes the effects of the timing model, WN, individual pulsar RN, 
if significant, and the GWB. In panel (c) P&(f) has been converted into units of 
characteristic strain, h.(f), for easy comparison with other GW detectors. 


highest frequencies. (Note that any DM variation model with 
the same frequency resolution would remove a similar 
amount of power.) 
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Figure 9. Stochastic GWB sensitivity curves. The solid blue curve shows the characteristic strain sensitivity as a function of GW frequency for the NANOGrav 15 yr 
data set and includes contributions from the WN, RN, and GWB self-noise. The gray line shows the recovered GWB (NG15gwb), while the dashed orange curve is the 
sensitivity curve only including the contributions from WN. The importance of including time-correlated noise can be seen in the overly optimistic orange dashed 
curve. The vertical lines show frequencies at one and two times the inverse of the data set time span. 


The next step in calculating a pulsar sensitivity curve is to 
calculate V -!, the noise-weighted transmission function. 


Ns 


YiG(GTCO)' !GT jj ef Git) (21) 


span ijj 


which can also be thought of as the timing-model-marginalized 
power spectral density, i.e., the Fourier domain equivalent of 
the timing-model-marginalized covariance matrix. In fact, 
N = Pr(f) is the power spectral density of the pulsar's 
residuals, and hence the power output of a single arm of our 
Galactic-scale GW detector. The power spectral density can be 
transformed to units of GW strain by taking into account the 
response function of TOAs to GWs. 

Figure 8(b) shows three separate curves in residual power. 


The dotted line shows the WN power spectral density, PP", 


where PAWN) = > JOf M giri), This only includes the 
span i 


WN described in Section 3.1 and shows what the WN power 
would look like without the effects of the transmission 
function. The orange dashed line shows the effects of the 
transmission function on the WN power. It is calculated by 
replacing C in Equation (21) with C™™. The solid blue line 
includes all of the modeled noise, modeling the same power as 
shown in the black solid line of Figure 3. 

The power spectral density in terms of pulsar residuals, Pp, 
can be converted into strain power spectral density, 


Pr(f) 


R 


by taking into account the pulsar response function to GWs. 
The sky-averaged response function useful for GWB char- 
acterization is simply 


S(f) = = 127°? PRA), (22) 


1 


R(f) = 1223/2" 


Q3) 
where a factor of three comes from the sky averaging and the 
Of) stems from the fact that we are working with timing 
residuals instead of Doppler shifts in the pulse frequency, 
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where hows = Af./f;. See Hazboun et al. (2019b) and Sesana 
et al. (2004) for more details about the pulsar residual response 
function. 

The strain power spectral density can then, in turn, be 
converted into units of characteristic strain, the units in which 
sensitivity curves are often plotted, via 
hf) = JfS(Q) = nf? SIRC). Figure 8(c) shows the 
strain power spectral density for J1909—3744 converted into 
units of characteristic strain. In these units there is a distinctive 
positive slope in the curve at higher frequencies that goes like 
~f 3. It is also evident that the GWB acts as a noise floor for an 
individual pulsar at some frequencies. It is important to note 
that the correlated Earth-term power in the GWB as measured 
at Earth (“Earth term") is the signal that we are searching for, 
but the power at the pulsar (“pulsar term") is uncorrelated 
across the pulsars (Cornish & Sesana 2013) and hence needs to 
be included in the noise budget and calculations of the 
sensitivity. 

The individual pulsar strain power spectral densities can be 
combined using the optimal statistic S/N in the frequency 
domain to construct a detection sensitivity curve for the 
NANOGrav PTA’° (Hazboun et al. 2019b), shown in Figure 9. 
The blue curve is the fiducial curve to compare the stochastic 
gravitational signals against since it contains the self-noise due 
to the background itself. The orange curve, which only includes 
WN, shows the importance of including the GWB and RN in 
calculating the sensitivity. The sensitivity curve in Figure 9 is 
calculated for any possible spectral shape of the stochastic 
GWB. In order to understand the detectability of that 
background, one needs to integrate the GWB power spectral 
density against the full "effective" power spectral density in the 
frequency domain. Following Thrane & Romano (2013), we 
determine the power-law-integrated sensitivity curve by 
calculating the amplitude of a power-law GWB we would 
detect with an S/N=5, as reported by the optimal statistic 


76 These curves are constructed using a statistic specific to a stochastic GWB. 
Curves for single, resolvable sources can also be constructed using a matched- 
filter statistic (Hazboun et al. 2019b). 
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Figure 10. Power-law-integrated sensitivity curve. The full stochastic 
sensitivity curve for PTAs combines the individual pulsar sensitivities using 
the denominator of the optimal statistic (Chamberlin et al. 2015; Hazboun 
et al. 2019b). As in Thrane & Romano (2013), one can integrate the sensitivity 
curve against a set of power-law stochastic backgrounds, here shown as the set 


of gray lines with a = -i to a = i and solve for the amplitude needed, 


ÁGwp. to meet a given threshold in S/N. The power-law-integrated sensitivity 
curve (solid black) is the maximum of the set of those lines at each frequency. 
The blue line shows the GWB sensitivity, as described in Figure 9, and the 
orange and green lines show the power-law GWBs detectable with S/N = 5 for 
two discrete values of o, described in the text. 


in NGl5gwb. In Figure 10 we show the detection sensitivity 
curve in blue and the power-law-integrated curve in black. 
Two particular spectral index values are highlighted: a = — 2/ 
3 (y= 13/3), the value theoretically expected for a population 
of binaries where the loss of energy to GWs dominates their 
evolution, and a= — 0.18 (y= 3.35), the median value for a 
obtained from the GWB search in NGl5gwb. 

The power-law-integrated calculation. gives us another 
important accounting of the noise budget for the NANOGrav 
detector. This analysis demonstrates that by folding together 
our knowledge of the noise characterization of the individual 
pulsars and the parameters for the common process, we can 
estimate the recovered amplitude, E) for the two values of 
the GWB spectral index used in NGl5gwb. 


6. Improvements and Future Outlook 
6.1. Near-future Wideband Noise Mitigation Strategies 


Over the past several decades, the bandwidths of the radio 
telescope receivers used for pulsar observations have steadily 
increased, with many receivers now capable of observing over 
more than one or even several octaves of bandwidth. Some 
examples of receivers used for PTA science are the Parkes— 
Murriyang  Ultra-Wideband Low Receiver (UWL, 
0.7-4.2 GHz; Hobbs et al. 2020b), the Effelsberg Ultra- 
Broadband Receiver (UBB; 1.3-6 GHz’*), and, for NANO- 
Grav, a 0.7-4.0GHz Ultra Wideband Receiver (UWBR; 
Bulatek & White 2020) being commissioned at the GBT, as 
well as the octave capability of CHIME (400-800 MHz; 
CHIME/Pulsar Collaboration et al 2021). Processing the data 
streams remains a “big data” problem in which raw voltage 
data are not desirable. 

For pulsar timing applications, data are first coherently 
dedispersed to remove DM delays in real time. Traditionally, 


77 [n order to make the full curve aesthetically complete, a wide range of 
spectral indices are used. 
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TOAs were then calculated over a set of discrete frequency 
channels. Comparing TOAs over wide frequency ranges then 
yielded deviations of DM from the fiducial value used in the 
dedispersion kernel. However, this traditional procedure 
ignored the pulse profile evolution apparent across wide 
frequency ranges, which is covariant with the DM. This has 
long been accounted for in NANOGrav timing models through 
fiting frequency-dependent (FD) pulse profile parameters to 
the channelized TOAs. The move toward even wider and 
continuously sampled bandwidths, however, favors a newer 
method of timing TOAs, by representing the profile as a 2D 
function that takes into account a fiducial DM, pulse profile 
frequency evolution, and scattering in a procedure known as 
"pulse portraiture" (Liu et al. 2014; Pennucci et al. 2014), 
which reduces covariances. Starting with the 12.5 yr data 
release, NANOGrav has been creating both narrowband and 
wideband data sets (Alam et al. 2021a, 2021b). With each new 
data release, we expect the wideband data set to become more 
and more important, as data from the newer wideband receivers 
will become a larger fraction of each existing data set. 
Additionally, the portraiture framework will be critical for 
probing frequency-dependent DM behavior, expected for some 
pulsars across very wide bandwidths (Cordes et al. 2016). 

Not all near-future noise mitigation strategies depend on new 
hardware. We highlight a signal processing technique known as 
cyclic spectroscopy (CS). Demorest (2011) and Walker et al. 
(2013) used CS as a deconvolution method to analyze data 
from PSR B19374-21, both resolving an intrinsic profile and 
revealing scintillation features with better resolution than that 
derived from filter bank techniques. Additional work since has 
shown how we can improve pulsar scintillation and TOA 
measurements (Archibald et al. 2014; Palliyaguru et al. 2015; 
Dolch et al. 2021; Turner et al. 2023), where CS is most 
effective for pulsars with high-profile S/N and/or high 
scattering timescales, as is the case for PSR B1937+21. While 
PSR B1937+21 is more highly scattered than most current 
PTA pulsars, future telescopes and upcoming receivers on 
existing telescopes will likely discover new, distant pulsars for 
which CS deconvolution is applicable. In addition, some 
currently known pulsars may become PTA quality. 

We expect significant S/N improvement both from the 
highly resolved pulse scintillation structure enabled by CS and 
from improved RFI mitigation due to frequency binning. This 
will be especially important for the heavy RFI contamination 
expected in the UWBR context. The computational limit, 
however, is that CS requires either saved baseband data or 
cyclic spectra written in real time, much like pulse profiles are 
currently written in real time after the application of coherent 
dedispersion. For the second (preferable) case, a real-time CS 
back end on the GBT’s UWBR is currently under development. 


6.2. Noise Mitigation for Different GW Sources 


PTA collaborations continually strive to understand and 
mitigate the noise in pulsar data sets. A deep understanding of 
the noise allows PTA collaborations to forecast the sensitivity 
of our detector to various GW signals and even tune it toward 
specific sources in some cases (Liu et al. 2023; Speri et al. 
2023). Figure 9 clearly demonstrates the strong effect of the 
stochastic GWB on the sensitivity of our detector. As PTAs 
enter the post-GWB detection era, the GWB will in fact be a 
foreground source of noise, relative to the many individual GW 
sources searched for by PTAs. While continued resolution of 
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the background (Pol et al. 2021) will allow us to mitigate up to 
half of the power in the GWB, half of the power in any given 
pulsar is uncorrelated across the pulsars, i.e., the pulsar term. 
Flagship GW analyses, e.g., for resolvable SMBHBs (Arzou- 
manian et al. 2023), already necessitate joint searches for a 
GWB and the deterministic signal to mitigate false positives. 
Searches for steeper spectral index stochastic backgrounds will 
require continued long-time-span observations (Kaiser et al. 
2022), while searches for resolvable binary signals will be 
served by better sensitivity at higher frequencies (Liu et al. 
2023). 

In recent years a number of new strategies for noise 
modeling of pulsar timing data sets have been developed (Lam 
et al. 2018a; Goncharov et al. 2021a; Chalumeau et al. 2022) to 
mitigate against higher-order chromatic propagation effects like 
scattering and telescope issues (i.e., “band noise" and “system 
noise"). The latter noise types are characterized using power- 
law RN models exclusively in a particular frequency band or 
for a specific pulsar back end. Largely these advances come in 
the form of Gaussian process models, built into the Bayesian 
noise analyses and also used in the full PTA GW searches. In 
Falxa et al. (2023) it was shown that including a large number 
of frequencies in the DM modeling was necessary to mitigate 
against spurious high-frequency single-source detections. 
These efforts reveal the need for and potential of these noise 
mitigation strategies for increasing sensitivity to high-fre- 
quency, deterministic GW sources. 


6.3. The Future of PTA Noise Mitigation 


New noise mitigation strategies for PTAs offer a twofold 
benefit for the low-frequency search for GWs. Mitigation 
techniques that involve new data will allow for better and more 
sensitive searches of GWs from sources in the higher end of 
our frequency band. Strategies that allow for mitigation of 
noise across the 15-- yr data set allow us to extract more 
sensitivity out of the data that we already possess. These 
techniques, as well as the in-depth characterization of PTAs as 
Galactic-scale detectors, will allow PTAs to reach their full 
potential as sensitive nanohertz GW observatories. 


7. Conclusion 


NANOGrav is a Galactic-scale GW detector made up of 
some of the most sensitive radio telescopes in the world and 
nature's most stable clocks. The astrophysical makeup of our 
detector requires careful noise analysis, starting from the 
underlying physical processes adding uncertainty to our 
measurements. It also requires that we use a well-motivated 
phenomenological model in order to incorporate all sources of 
noise, whether explicitly accounted for or not, in order to be as 
conservative as possible in our characterization. 

We have carried out a full characterization of the noise in our 
nanohertz GW detector, including a full accounting of RN in 
the pulsars and new analyses to corroborate the use of power- 
law RN models. These noise analyses allow us to construct 
detection sensitivity curves for the stochastic GWB that are in 
agreement with the results of NGl5gwb. We find that our 
phenomenological model faithfully encapsulates the noise in 
our detector. 
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graphy. J.S.H., N.La., J.S., and X.S carried out the noise 
analyses. J.S.H carried out the sensitivity analysis. 
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Appendix A 
The Noise Covariance Matrix 


A.1. Marginalizing the Timing Models 


In Equation (9), we identify the full PTA marginalized 
timing model as 


P(6t\a, E) = ILC a, E)P(e)d(e) (Al) 
1 T p-1 
exp( — 1 C&t — Fa) D` (ôt — Fa)) 
= ; (A2) 
ydet(27D) 
where 
D = C + MXM’. (A3) 
Using the Woodbury identity, D^! can be written as 
D-!=(C + MXM’)! 
-Cl-CM (QX-!-- MT'C-MyMTC- (A4) 


One can now set the range on the priors for e to be infinite by 
setting X = diag( oo); hence, X ! — 0 and the expression in 
Equation (A4) is further simplified to 


D-ciR-c!-cC-MM(M'C-My!MTC- (AS) 
which has completely removed any dependence on the timing 
model parameters from the likelihood, though dependence on 
the timing model is preserved in the design matrix, M. Note 
that we have defined R = 1 — M(MTC-M)-1M'C-|, where 
R acts to project the effects of the timing model out of the 
covariance matrix. 

We use a similar marginalization to find the final form for a 
single-pulsar likelihood (Equation (13)): 


P&E) = f P(6tla, E)P(a)d (a) 


B exp(— zo N-ter) 


(A6) 
A det(2xN) 
where 
NzD-4 FoF’. (AT) 


The Woodbury lemma can again be used to get the following 
form for the inverse: 


N-!« (D + FoF’)! 
= D! — DIF (y7! + FTD-IFy FTD! (A8) 
The marginalizations discussed above and used to define the 
noise covariance matrix in Section 4.1 have been developed 
precisely for efficient calculation since D™' and F are constant 
when the WN is fixed, and the only matrix inversion that needs 
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to be calculated is the one in the parenthetical of Equation (A8), 
a matrix that is only (ZNgeq X 2Nfreq): 


A.2. Expanded Signal Covariance Matrix 


As described in Section 4, we define the noise model of a 
particular observing epoch as the matrix containing information 
about the frequency dependence of template-fitting uncertain- 
ties asyn, pulse jitter oj, and diffractive interstellar scattering 
Opiss (see Sections 3 and 2, respectively): 


2 2 2 2 
Os, Q^) + OF + Opiss(M1, vı) 


2 2 
epoch — 0; + Opiss(Y2, vi) 


2 2 
og + C piss v; Vi) 


Compare Equation (A9) to the phenomenological covariance 
matrix presented in Equation (14). The diagonal elements 
describe the autocovariance at each observing frequency and 
hence include cs, whereas the off-diagonal elements describe 
noise correlated across frequencies, here jitter and diffractive 
interstellar scattering. Using this covariance matrix for each 
observing epoch, we can construct a full noise model following 
the structure of Equation (15): 


2 
0; + Opiss(i v2) 


2 2 2 
og /n(¥2) + 05 + Gbiss(2. V2) + 


2 2 
og +t C piss YN, V2) 
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epoch ( A) 0 NA 0 
C= 0 CePoch (t) -— 0 
0 0 Cepoch ( ty) 


Hn C observatory E (C 65M zd CsPin (|t == Vv), (A10) 


where CP?” is modeled as a power-law spectral process and 
CoPservatory ang C°>™ are derived directly from observation and 


2 2 
07 + Opiss (i, VN) 


oF EN OG orss (V2; VN) (A9) 


2 2 2 
OS NN) + 07 + Opiss(/v. VN) 


encapsulate noise from the radio telescope and variations in the 
ISM, respectively. 


Appendix B 
Definitions 


In Table 3 we include an exhaustive list of the notation used 
throughout this work. A carat () denotes an estimator, and we 
have excluded certain subscripts on variables listed here 
(e.g., GWB). 
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Table 3 
Symbols and Acronyms Used 


Symbols Definition Typical Units 
a Fourier amplitudes vector 
A Spectral amplitude Strain 
(unitless) 
B Observing bandwidth MHz 
B Bayes factor 
C Generic covariance matrix 
Cree Covariance matrix of white noise on given 
epoch 
cobservatory ^ Covariance matrix for noise at the observatory 
CoM Covariance matrix for DM noise 
crn Covariance matrix for spin noise 
CEN Red-noise covariance matrix 
df Fourier frequency differential for power Hz 
spectral density integrand 
D D =C + MXM" 
DM Dispersion measure pe cm~ 
E Set of white-noise parameters (F, Q, J) 
f Spectral /Fourier frequency Hz 
S Nyquist frequency 
te Pulsar spin frequency Hz 
te Frequency of 1 yr! yr! 
F Fourier design matrix 
G [G matrix] 
F EFAC, multiplicative error factor 
h Strain 
he Characteristic strain 
i,j Time-of-arrival indices 
"d ECORR, error correlated across radio us 
frequencies 
M Design matrix 
n White noise (vector) 
no Electron number density at a distance of 1 au cm? 
N N=D + FoF" 
Nireq Number of Fourier frequencies 
N, Number of pulses 
Nroa Number of times of arrival 
N Power spectral density of pulsar resi- 
duals, N = Pr(f) 
Note that M~! is the noise-weighted trans- 
mission function 
P(f) Power spectral density Hz! 
PRO) Power spectral density of residuals sHz! 
PWN White-noise power spectral density of the 
residuals 
P, Pulsar spin period ms 
P Probability density function 
p Model parameters (vector) 
Po Best-fit model parameters (vector) 
Q EQUAD, error added in quadrature us 
r Radio frequency ratio 
RO) Sky-averaged response function for the GWB 
S Strain power spectral density 
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S/N 


Table 3 
(Continued) 

Definition 

Signal-to-noise ratio 

Time of arrival 

Times of arrival (vector) 

Timing model transmission function 

Nondispersive chromatic delays 

Scattering delay timescale 

Separation time between multifrequency 
observations 

Single observation length 

Total time span of data set 

Block-diagonal matrix with values of 1 for 
TOAs from 

the same observation epoch and 0 otherwise 

Pulse width 

Covariance matrix of timing model 
parameters 

Spectral index of GWB in units of character- 
istic strain 

Ecliptic latitude 

Spectral index 

Kronecker delta 

Timing residuals (vector) 

Deterministic non-GW model times of arrival 
(vector) 

Model/predicted times of arrival (vector) 

Perturbation on the infinite-frequency arri- 
val time 

Diffractive scintillation timescale 

Change in pulsar spin frequency 

Diffractive scintillation bandwidth 

Perturbations to the model parameters (vector) 

Frequency-dependent timing errors 

Fractional gain error 

Radio frequency 

Degree of circular polarization 

Diffractive interstellar scintillation noise 

Jitter noise 

Polarization miscalibration uncertainty 

Signal-to-noise-ratio-dependent template-fit- 
ting uncertainty 

Spin noise 

Infinite-frequency arrival time uncertainty 

Sun—pulsar separation angle 

Variance of RN coefficients (see 
Equation (11)) 
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Typical Units 


Ls 


min 
yrs 


deg 


Ls 


MHz, GHz 


Ls 
Ls 


deg 


Noise budget figures for all of the pulsars used in Agazie et al. 
(2023b) are provided in Figures 11-73. See Section 5.2.1, and the 
caption of Figure 3 for details. 
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Figure 11. The excess timing residual delay as a function of frequency for PSR B1855+09. See Figure 3 for details. 
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Figure 12. The excess timing residual delay as a function of frequency for PSR B1953+29. See Figure 3 for details. 
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Figure 13. The excess timing residual delay as a function of frequency for PSR J0023+0923. See Figure 3 for details. 
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Figure 14. The excess timing residual delay as a function of frequency for PSR J0030--0451. See Figure 3 for details. 
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Figure 15. The excess timing residual delay as a function of frequency for PSR J0340--4130. See Figure 3 for details. 
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Figure 16. The excess timing residual delay as a function of frequency for PSR J0406--3039. See Figure 3 for details. 
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NG15yr Noise Budget: J0437—4715 (Tipan=4.8 yrs) 
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Figure 17. The excess timing residual delay as a function of frequency for PSR J0437—4715. See Figure 3 for details. 
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Figure 18. The excess timing residual delay as a function of frequency for PSR J0557--1551. See Figure 3 for details. 
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Figure 19. The excess timing residual delay as a function of frequency for PSR J0605+3757. See Figure 3 for details. 
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Figure 20. The excess timing residual delay as a function of frequency for PSR J0610—2100. See Figure 3 for details. 
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Figure 21. The excess timing residual delay as a function of frequency for PSR J0613—0200. See Figure 3 for details. 
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Figure 22. The excess timing residual delay as a function of frequency for PSR J0636--5128. See Figure 3 for details. 
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Figure 23. The excess timing residual delay as a function of frequency for PSR J0645--5158. See Figure 3 for details. 
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Figure 24. The excess timing residual delay as a function of frequency for PSR J0709--0458. See Figure 3 for details. 
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Figure 25. The excess timing residual delay as a function of frequency for PSR J0740--6620. See Figure 3 for details. 
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Figure 26. The excess timing residual delay as a function of frequency for PSR J0931— 1902. See Figure 3 for details. 
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Figure 27. The excess timing residual delay as a function of frequency for PSR J1012--5307. See Figure 3 for details. 


NG15yr Noise Budget: J1012—4235 (Tipan=3.4 yrs) 


9€ Excess Noise Posteriors **** [Individual Red -——- Total 
9 Free Spectral Posteriors White Noise == = CURN 


A ——— — Lec 


i. —— =~ — — {— f -+t 


Frequency (Hz) 


Figure 28. The excess timing residual delay as a function of frequency for PSR J1012—4235. See Figure 3 for details. 
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Figure 29. The excess timing residual delay as a function of frequency for PSR J1022-- 1001. See Figure 3 for details. 
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Figure 30. The excess timing residual delay as a function of frequency for PSR J1024—0719. See Figure 3 for details. 
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Figure 31. The excess timing residual delay as a function of frequency for PSR J1125--7819. See Figure 3 for details. 
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Figure 32. The excess timing residual delay as a function of frequency for PSR J1312+0051. See Figure 3 for details. 
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Figure 33. The excess timing residual delay as a function of frequency for PSR J1453--1902. See Figure 3 for details. 
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Figure 34. The excess timing residual delay as a function of frequency for PSR J1455—3330. See Figure 3 for details. 
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Figure 35. The excess timing residual delay as a function of frequency for PSR J1600—3053. See Figure 3 for details. 
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Figure 36. The excess timing residual delay as a function of frequency for PSR J1614—2230. See Figure 3 for details. 
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Figure 37. The excess timing residual delay as a function of frequency for PSR J1630--3734. See Figure 3 for details. 
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Figure 38. The excess timing residual delay as a function of frequency for PSR J1640--2224. See Figure 3 for details. 
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Figure 39. The excess timing residual delay as a function of frequency for PSR J1643— 1224. See Figure 3 for details. 
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Figure 40. The excess timing residual delay as a function of frequency for PSR J1705— 1903. See Figure 3 for details. 
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Figure 41. The excess timing residual delay as a function of frequency for PSR J1713+0747. See Figure 3 for details. 
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Figure 42. The excess timing residual delay as a function of frequency for PSR J1719—1438. See Figure 3 for details. 
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Figure 43. The excess timing residual delay as a function of frequency for PSR J1730—2304. See Figure 3 for details. 
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Figure 44. The excess timing residual delay as a function of frequency for PSR J1738--0333. See Figure 3 for details. 
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NG15yr Noise Budget: J17414-1351 (Tspan=11.0 yrs) 
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Figure 45. The excess timing residual delay as a function of frequency for PSR J1741--1351. See Figure 3 for details. 
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Figure 46. The excess timing residual delay as a function of frequency for PSR J1744— 1134. See Figure 3 for details. 
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Figure 47. The excess timing residual delay as a function of frequency for PSR J1745--1017. See Figure 3 for details. 
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Figure 48. The excess timing residual delay as a function of frequency for PSR J1747—4036. See Figure 3 for details. 
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NG15yr Noise Budget: J1751—2857 (1,54, —3.5 yrs) 
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Figure 49. The excess timing residual delay as a function of frequency for PSR J1751—2857. See Figure 3 for details. 
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Figure 50. The excess timing residual delay as a function of frequency for PSR J1802—2124. See Figure 3 for details. 
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Figure 51. The excess timing residual delay as a function of frequency for PSR J1811—2405. See Figure 3 for details. 
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Figure 52. The excess timing residual delay as a function of frequency for PSR J1832—0836. See Figure 3 for details. 
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Figure 53. The excess timing residual delay as a function of frequency for PSR J1843—1113. See Figure 3 for details. 
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Figure 54. The excess timing residual delay as a function of frequency for PSR J1853--1303. See Figure 3 for details. 
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Figure 55. The excess timing residual delay as a function of frequency for PSR J1910--1256. See Figure 3 for details. 
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Figure 56. The excess timing residual delay as a function of frequency for PSR J1911+1347. See Figure 3 for details. 
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Figure 57. The excess timing residual delay as a function of frequency for PSR J1918—0642. See Figure 3 for details. 
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Figure 58. The excess timing residual delay as a function of frequency for PSR J1923+2515. See Figure 3 for details. 
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Figure 59. The excess timing residual delay as a function of frequency for PSR J1944--0907. See Figure 3 for details. 
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Figure 60. The excess timing residual delay as a function of frequency for PSR J1946--3417. See Figure 3 for details. 
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Figure 61. The excess timing residual delay as a function of frequency for PSR J2010— 1323. See Figure 3 for details. 
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Figure 62. The excess timing residual delay as a function of frequency for PSR J2017+0603. See Figure 3 for details. 
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Figure 63. The excess timing residual delay as a function of frequency for PSR J2033--1734. See Figure 3 for details. 
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Figure 64. The excess timing residual delay as a function of frequency for PSR J2043+1711. See Figure 3 for details. 
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Figure 65. The excess timing residual delay as a function of frequency for PSR J2124— 3358. See Figure 3 for details. 
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Figure 66. The excess timing residual delay as a function of frequency for PSR J2145—0750. See Figure 3 for details. 
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Figure 67. The excess timing residual delay as a function of frequency for PSR J2214+3000. See Figure 3 for details. 
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Figure 68. The excess timing residual delay as a function of frequency for PSR J2229--2643. See Figure 3 for details. 
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Figure 69. The excess timing residual delay as a function of frequency for PSR J2234--0611. See Figure 3 for details. 
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Figure 70. The excess timing residual delay as a function of frequency for PSR J2234+0944. See Figure 3 for details. 
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Figure 71. The excess timing residual delay as a function of frequency for PSR J2302+4442. See Figure 3 for details. 
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Figure 72. The excess timing residual delay as a function of frequency for PSR J2317+1439. See Figure 3 for details. 
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Figure 73. The excess timing residual delay as a function of frequency for PSR J2322+2057. See Figure 3 for details. 
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